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ABSTRACT 

Supermassive black hole (BH) mergers produce powerful gravitational wave (GW) emission. 
Asymmetry in this emission imparts a recoil kick to the merged BH, which can eject the 
BH from its host galaxy altogether. Recoiling BHs could be observed as offset active galac¬ 
tic nuclei (AGN). Several candidates have been identified, but systematic searches have been 
hampered by large uncertainties regarding their observability. By extracting merging BHs 
and host galaxy properties from the Illustris cosmological simulations, we have developed 
a comprehensive model for recoiling AGN. Here, for the first time, we model the effects of 
BH spin alignment and recoil dynamics based on the gas-richness of host galaxies. We pre¬ 
dict that if BH spins are not highly aligned, seeing-limited observations could resolve offset 
AGN, making them promising targets for all-sky surveys. For randomly-oriented spins, ^10 
spatially-offset AGN may be detectable in HST-COSMOS, and > 10^ could be found with 
Pan-STARRS, LSST, Euclid, and WFIRST. Nearly a thousand velocity-offset AGN are pre¬ 
dicted within the SDSS footprint; the rarity of large broad-line offsets among SDSS quasars is 
likely due in part to selection effects but suggests that spin alignment plays a role in suppress¬ 
ing recoils. Nonetheless, in our most physically motivated model where alignment occurs only 
in gas-rich mergers, hundreds of offset AGN should be found in all-sky surveys. Our findings 
strongly motivate a dedicated search for recoiling AGN. 

Key words: accretion, accretion discs - black hole physics - gravitational waves - hydrody¬ 
namics - galaxies: active - galaxies: interactions 


1 INTRODUCTION 

Mergers between galaxies are important drivers of evolution. They 
grow stellar bulges, evolve galaxies along the Hubble sequence, 
and can trigger both star formation and accretion onto central su¬ 
permassive black holes (BHs), which in turn influence the host via 


feedback processes (e.g., Sanders et al.|1988| Barnes & Hemquist 

1991||1996[ 

Mihos & Hemquist||1996l|Wyithe & Loeb||2003 |Di 

Matteo et al. 

2005||Hopkins et al.|2006||2008a|l. While the majority 


of active galactic nuclei (AGN) in the local Universe do not appear 


to be triggered by me rgers dCisternas et al.|20ir||Schawinski et al.| 
|2011[ |Kocevski et al.||2012| l, a clear correlation between merging 
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activity and AGN fueling has been observed (e.g., |Koss et al.|2010[ 
[Ellison et al.|201 l|[Satyapal et al.|2014| l. There is also evidence that 
the majority of the highest-luminosity quasars are triggered by ma¬ 
jor mergers (e.g., jUrmtia et al.|2008||Treister et al.|2012||Hopkins| 
jet al.|2014| l. 

In addition to fueling rapid growth, a major merger between 
comparable-mass galaxies, each containing a central BH, will in¬ 
evitably lead to the formation of a BH binary. This binary may 
eventually merge, releasing vast amounts of energy as gravitational 
waves (GWs). Such events could be detected directly in the com¬ 
ing years via pulsar timing arrays (PTAs, e.g., |Sesana et al.|2008[ 
|Sesana|2013j ), or with a future space-based GW interferometer such 
as eLISA (e.g., |Amaro-Seoane et al.|20T^ . Additionally, any asym¬ 
metry in the merging BH system - unequal masses or spins - will 
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2 L. Blecha et al. 


produce asymmetric GW emission and a net momentum flux, caus¬ 
ing the merged BH to “recoil” in the opposite direction at the mo¬ 
ment of merger ( |Peres|19^|Bekenstein|1973] l. Numerical relativ¬ 
ity simulations in recent years have demonstrated that these kicks 
may be very large, up to 5000 km s“^, implying that some BHs 
could be ejected from galaxies entirely ( [Campanelli et aL||2007[ 
[Lousto & Zlochower|2011| l. 

These “superkicks” of more than a thousand km s“^ should 
be rare, but smaller recoils may still displace BHs from galactic 
nuclei for substantial periods of time. Recoil events may have been 
more frequent in the early Universe, when galactic escape speeds 
were lower and merger rates were higher. This places extra con¬ 
straints on models of BH seed formation that must produce massive 
quasars by z ~ 7 (e.g., |Volonteri & Rees|2006[|Li et al.|2007[|Si-| 
[jacki et al.|2009| l. At lower redshift, BH ejections must be relatively 
rare in massive ellipticals, at least, owing to the observed ubiquity 
of central BHs in such galaxies for which mass measurements have 
been obtained. The sample of galaxies with direct mass measure¬ 
ments is still modest, however, and the lower-mass range is largely 
unexplored. A large sample of dynamical mass measurements in 
massive galaxies, such as brightest cluster galaxies (BCGs), could 
yield strong evidence for a past recoil event if a galaxy with a 
“missing” or very undermassive BH were found ( [Gerosa & Sesana| 
|2015| ). Even if the BH is not ejected entirely, a marginally-bound 
BH may oscillate in the galactic potential for many Gyr before re¬ 
turning to the center (e.g.,[Madau & Quataert|2004| Gualandris &| 


|Merritt|2008[|Komossa & Merritt|2008[[Blecha & Loeb|2008| l. 

These GW recoil events have implications for BH-galaxy co¬ 
evolution. They increase the intrinsic scatter to the observed BH- 
bulge relations ( |Volonteri|2007|l, even wh en the BH is not ejected 
entirely from the host ( Blecha et al.|2011||Sijacki et al.|2011| l. The 
displacement of the BH from the galactic center also reduces the 
amount of AGN feedback imparted to the nuclear region, allowing 
a denser stellar cusp to form and delaying the galaxy’s transition 
from the blue to the red sequence ( [Blecha et al.|2011| ). 

A major uncertainty in predicting the effects of GW recoil is 
that the kick velocity distribution depends sensitively on the pre¬ 
merger BH spin vectors, which are not known. In particular, if the 
BH binary is embedded in a circumbinary gas disk that drives the 
binary evolution, torques from the disk may efficiently align the BH 
spins prior to a merger (e.g., |Bogdanovic et al.||2007[|Dotti et al.j 
|2010[ [Miller & Krolik||2013j ). The maximum recoil kick for per¬ 

fectly aligned spins is < 200 km s“^, even for maximally spinning 
BHs, and kicks > 600 km s“^ do not occur if the spins are aligned 
to within a few degrees of the orbital angular momentum. However, 
even in gas-rich systems where a circumnuclear disk forms, spin 
alignment may not always be efficient ( [Lodato & Gerosa[[2013[ l. 
Also, if the accretion flow is chaotic, characterized by low-mass 
gas streams (< 1% Mbh) infalling from random directions, the 
BH spins may not be efficiently aligned prior to merger. In cer¬ 
tain cases, circumbinary disks could even produce stable counter- 
alignment if the disk angular momentum is low and the initial BH 
misalignment is large ( [Nixon[20'T^ , and retrograde disks may drive 
the BHs to coalescence before alignment can occur ( [Schnittman &[ 
[Krolik[20'T5] l. Such configurations may occur if outer regions of the 
disk become unstable to star formation, or if the nuclear gas inflow 
is intermittent (e.g., [King & Pringle[2006[[Lodato et al.[2009[[King[ 
[& Nixon|2013| l. Furthermore, stellar dynamics may be the domi¬ 
nant driver of BH binary inspiral in some systems, especially in tri- 
axial, gas-poor mergers (e.g.,[Berczik et al.[2006[[Khan et al.[201 1| 


[Holley-Bockelmann & Khan|2015[ l. These stellar interactions typi¬ 

cally will not cause any preferential alignment of BH spins ( [Merritt[ 


[& Vasiliev|2012[ |. In all cases, even in the final stages of BH inspi¬ 
ral after the BHs have dynamically decoupled from their surround¬ 
ings, their spins will undergo further evolution via relativistic pre¬ 
cession. This precession can substantially alter the final BH spins 
prior to merger, driving them toward alignment in some cases and 
increasing misalignment in others, thereby decreasing or increas¬ 
ing the resultant recoil velocity depending on the mass ratio and 
initial spins of the BH binary ( [Berti et al.[2012[[Kesden et al.[2015[ 
[Gerosa et al.|2015[ l. 

In addition to its possible role in aligning BH spins, gas 
in merging galaxies may strongly influence the recoil trajectories 
themselves. The inflow of cold gas to central regions during a 
gas-rich major merger can dramatically increase the central escape 
speed of the galaxy for a period of time, and gas drag can also 
help retain BHs in galactic nuclei ([Blecha et al.|2011|[Sijacki et al.[ 
|20TT1 ). Thus, even if spin alignment via gas dynamics is less ef¬ 
ficient than numerical results suggest, large BH displacements in 
gas-rich systems may be rare. 

In contrast, recoil events in gas-poor mergers can be much 
longer lived, particularly if the inner stellar profile has a shallow 
density core. Cuspy stellar density profiles may be transformed 
into cores (and cores may be enlarged) by BH binary inspiral (e.g., 
[Quinlan|1996|[Yu|2002[[Milosavljevic & Merritt|2001[ | and by the 
impulsive removal of the BH and its bound cusp of stars ([Merritt[ 
[et al.|2004| l. 


A kicked BH that is actively accreting will carry along its in¬ 
ner accretion disk and broad emission line (BL) region, and may 
also sweep up gas from its surroundings, such that it could be ob¬ 
served as an AGN that is spatially- or kinematically-offset from its 

host galax}^(e.g., Madau & Quataert|2004 Loeb|2007 Blecha & 

- - 


[Loeb|200^ [Volonteri & Madau|2008[ Fujita|2009[ ). In the case of 

kinematic offsets, the line-of-sight (LOS) velocity of the recoiling 
BH could be detected in the object’s spectrum as an offset between 
the BLs and the narrow emission lines (NLs), which will remain 
bound to the host galaxy (e.g., [Merritt et al.[2006[ l. An offset could 
also be detectable between the AGN BLs and the rest-frame red- 
shift of the host galaxy itself. Observable offset AGN lifetimes can 
be up to tens of Myr for both spatial and kinematic offsets, for a 
wide range in kick speeds ( [Blecha et al.[2011[[Guedes et al.[201l] |. 
Prior to the direct detection of GWs, recoiling AGN could offer the 
best evidence of BH mergers. 

In the last few years, a handful of candidate (spatially- or 


velocity-) offset AGN have been identified 

([Komossa et al. 

2008 

Shields et al.|2009b|[Comerford et al.|2009 

[Robinson et al. 

2010 

Batcheldor et al.|20101[Civano et al.|2010| 

Steinhardt et al. 

2012 

Koss et al.|2014|[Lena et al.|2014|). None have yet been confirmed 


favored for some ( [Shields et al.[[2009¥| [Decarli et al.[[2014] ), but 
two candidates in particular appear promising. CID-42 is a highly 
disturbed galaxy with an AGN that appears both spatially and kine¬ 
matically offset from the host nucleus, though the possibility that 
the galaxy is in a pre-merger phase, with a second, quiescent BH in 
the nucleus, is difficult to rule out ( Comerford et al.|2()0^ [Civano[ 


let al.|2010|[20T2l|Blecha et al.|2013[|Novak et al.|2015t . SDSSl 133 


is a bright point source offset from a nearby dwarf galaxy; it has 
no evidence for an extended emission component down to a scale 
of 12 pc ( [Koss et al.[2014[ ) and is observed over a 65-year baseline. 


^ Following convention, we refer to any actively accreting BH as an “AGN” 
but note that in the event of a GW recoil, active BHs will generally not be 
coincident with the galactic nucleus. 
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Detection of recoiling BHs 3 


An alternative scenario in which the object experienced an extreme, 
50-year luminous blue variable type outburst followed by an unusu¬ 
ally long-lived Type Iln supernova cannot be excluded, but further 
monitoring of the source should be conclusive in the near future. 

These largely serendipitous discoveries motivate a systematic 
search for offset AGN. To date, a few searches for velocity-offset 
quasars have been carried out ( Bonning et al.|[2007l |Tsalmantza| 


|et al.|20lT||Eracleous et al.|2012| l. |Bonning et al" ( 20()7| l searched 

for consistent, symmetric velocity shifts in broad H/3 and Mg II 
lines in SDSS DR5 quasar spectra and foun d a null result for 
Au > 800 km s“ 


Tsalmantza et al. 


<201 it searched for offset 


BLs among SDSS DRV quasars, focusing on binary BH candidates, 
and identified 32 quasars with significant BL shifts. |Eracleous et ^ 
^20\2\ similarly examined the SDSS DRV quasars for broad H/3 
lines offset by Av > 1000 km s“^ and found 88 such objects. In 
these samples, some objects have BL shifts > 5000 km s“^, the 
theoretical maximum recoil velocity, and some are best explained 
as double-peaked emitters. Offset BLs may also be produced by 
high-velocity outflows or the orbital motion of a single active BH in 
a sub-parsec binary pair. These results indicate that recoiling AGN 
with observable velocity offsets are rare among SDSS quasars, if 
present at all. 

Spatially-offset AGN remain largely unexplored. [Lena et aT] 
( |2014| ) have conducted the only systematic search for spatial offsets 
thus far, focusing on very small-scale offsets (< 10 pc) in nearby 
giant core ellipticals. They find displacements between the AGN 
and the galactic photocenter in 6 of 14 ellipticals, including M8V 
where a displacement was previously noted by [Batche ldor et al.| 
( [20 To 1 ). However, the alignment of the displacement axis with a 
radio jet in four of the AGN argues against the recoil scenario for 
these objects. No searches have yet been carried out for AGN with 
larger spatial offsets. 

Designing a targeted search for offset AGN requires better 
theoretical understanding of when and where significant displace¬ 
ments are most likely to occur. Using semi-analytic models, and 
assuming randomly-oriented BH spins, |Volonteri & Madau| ( [2008| 
hereafter VM08) predict rates of observable spatially-offset AGN; 
their results are of particular relevance to this work and are com¬ 
pared with our findings in Section |4.2| The effect of pre-merger 
spin alignment remains a major open question, and the infiuence of 
gas on recoil dynamics has been studied only for isolated galaxy 
mergers, rather than in a cosmological framework. The host galaxy 
properties of offset AGN are also unexplored, with the exception 
of a recent study ( |Gerosa & Sesana|2015| l, which concludes that a 
larger sample of BH mass measurements in brightest cluster galax¬ 
ies (BCGs) could find evidence for past superkicks. 

Here, we address these open questions by constructing a 
model for the observability of spatially- and velocity-offset AGN, 
utilizing data from the state-of-the-art cosmological hydrodynamic 
simulation Illustris (e.g., |Vogelsberger et al.|2014a|b[ [Nelson et al.| 
|2015[[Genel et al.|2014| |. The Illustris simulation has been shown 

to successfully reproduce many observed properties of galaxies and 
their BHs, including the stellar and BH mass functions, cosmic star 
formation rate density, galaxy merger rate, baryonic Tully-Fisher 
relation, and quasar luminosity function ( [Vogelsberger et al.|2014b[ 
|Genel et al.|2014[[Sijacki et al.|2015| l. We examine the offset life¬ 
times and space density of observable offset AGN that could be 
detected in various surveys, and we make the first predictions re¬ 
garding the host galaxy properties of offset AGN. We consider 
a range of possible pre-merger BH spin distributions in order to 
study the effect of spin alignment on offset AGN observability, and 
we discuss how detections of recoils might be used to constrain 


the degree of alignment. These predictions will be testable in ex¬ 
isting and future surveys with the Hubble Space Telescope (HST), 
Chandra, the Panoramic Survey Telescope & Rapid Response Sys¬ 
tem (Pan-STARRS), the James Webb Space Telescope (JWST), the 
Large Synoptic Survey Telescope (LSST), Euclid, and the Wide- 
Field Infrared Survey (WFIRST). 

The outline of this paper is as follows. The cosmological sim¬ 
ulations are described in Section |2.1| BH spin models and recoil 
velocity distributions are detailed in Sec tion |2.2| a nd the model 
for recoiling AGN is outlined in Sections |2.3| - |2.5| Section de - 
scribes our results concerning the BH merger rate (Section [TT] ), 
recoil trajectories (Section [3^ , recoiling AGN properties (Section 
|3.3| ), spatial and velocity offsets (Sectionoffset AGN observ¬ 
ability (Section [33] ), and host galaxy properties (Section [33] ). The 
dependence of our results on model assumptions is discussed in 
Section [UVj We discuss the prospects for identification and follow¬ 
up of promising recoil candidates in Section jUT] and compare with 
previous work in Section [43] Our conclusions are summarized in 
Section [5] 


2 METHODOLOGY 

Our basic procedure for constructing a recoiling AGN model is as 
follows: 1) extract the properties of merging BH pairs and their 
host galaxies from the Illustris cosmological simulations, 2) assign 
pre-merger spins to the BHs and calculate the resulting recoil kick 
velocity, according to an assumed spin distribution, 3) construct an 
analytic potential model for the galaxy merger remnant at the time 
of BH merger, based on the properties of the progenitor galaxies, 4) 
integrate the trajectory of the recoiling BH in this potential, includ¬ 
ing stellar dynamical friction, and 5) calculate the AGN luminosity 
of the recoiling BH at each timestep and determine its observability 
as an offset AGN. Here we describe this procedure in detail. 


2.1 Cosmological simulations 

We use the data from the Illustris cosmological simulation project 
(e.g., [Vogelsberger et al.||2014b[ [Genel et3T][2014| > as the basis 
for our models. The Illustris simulations were conducted with the 
moving-mesh hydrodynamics code AREPO ( [Springel||2010j ). The 
highest-resolution simulation (“Illustris-1”) spans a comoving vol¬ 
ume of (106.5 Mpc)^ and has 2 x 1820^ resolution elements, 
with a dark matter (DM) mass resolution of 6.26 x 10® M© and 
a typical gas cell mass of 1.26 x 10® M©. The “Hlustris-2” and 
“niustris-3” simulations have the same volume with 2 x 910® and 
2 X 455® resolution elements, respectively. The following cosmo¬ 
logical parameters, which are consistent with the Wilkinson Mi¬ 
crowave Anisotropy Probe (WMAP)-9 measurements, are used in 
the simulation and assumed throughout this paper: Um = 0.2726, 
Qa = 0.7274, Qb = 0.0456, erg = 0.809, and Ho = lOOh km 
s“^ Mpc“^ with h = 0.704. Sub-resolution physical models are 
included for primordial and metal-line cooling, star formation and 
stellar feedback { [Springel & Hernquist|[2003| ), gas recycling and 
chemical enrichment, and BH seeding, accretion, and feedback. 
These prescriptions are described in detail in [Vogelsberger et al.] 
( [2013j ) and [Torrey et al.| ( [2014j ). Here, we outline the relevant as¬ 
pects of the prescriptions for BH formation, growth, and feedback 
in the Illustris simulation. These models successfully reproduce BH 
and AGN populations in good agreement with observational con¬ 
straints, including the BH mass function at z = 0 and the quasar 
luminosity function ( [Sijacki et al.|2015| l. 
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4 L. Blecha et al. 


Throughout the simulation, an on-the-fly Friends-of-Friends 
(FOF) algorithm is used to identify bound subhalos, and a BH parti¬ 
cle with a seed mass of M© is placed in each subhalo more 

massive than 7.1 x 10^° Mq that does not already contain a BH. 
Gas accretion onto BHs is parametrized in terms of a Bondi-like 
formula, capped at the Eddington limit. A repositioning scheme is 
implemented to center BHs on the gravitational potential minimum 
of their hosts, which prevents spurious BH wandering from numer¬ 
ical two-body scattering. As a result, BH velocities are ill-defined, 
and the velocity of the BH relative to the surrounding gas is not 
taken into account in calculating the accretion rate. Additionally, 
a pressure criterion is introduced to reduce the BH accretion rate 
if the gas pressure near the BH is insufficient to compress the gas 
against thermal feedback, which prevents the formation of unphys- 
ically large, hot gas bubbles around accreting BHs. 

In addition to gas accretion, BHs are allowed to grow via 
mergers. Two BHs are assumed to have merged when their sepa¬ 
ration falls below the gravitational softening length. Owing to the 
effect of the repositioning scheme on BH velocities, no criterion is 
imposed on the relative BH velocity in order for them to merge. 
Furthermore, in some instances, when a satellite halo containing a 
BH is determined to have merged with the central halo, the reposi¬ 
tioning scheme may cause the low-mass satellite BH to merge with 
the central BH on an artificially short timescale. We find that this 
effect is largely confined to BHs with mass < 10® M©, and we 
accordingly exclude from our analysis all mergers in which either 
BH mass is below this value. 

Three distinct modes of BH feedback are included in the sim¬ 
ulation. Thermal or “quasar mode” feedback is introduced by cou¬ 
pling 5% of the AGN luminosity to surrounding gas particles as 
thermal energy ( [Springel et al.|2005] l. When the accretion rate falls 
below 0.05 of the Eddington rate, AGN feedback is assumed to be 
dominated by a radiatively inefficient “radio mode,” in which hot, 
buoyant bubbles with a radius of 100 kpc are injected into the am¬ 
bient gas { [Sijacki et al.||2007| ). Additionally, “radiative” feedback 
is included by assuming a fixed AGN spectral energy distribution 
(SED) and modifying the net gas cooling rate in the presence of 
strong ionizing radiation. This radiative feedback is most relevant 
for BHs accreting near the Eddington limit. 

In addition to this constraint on the BH mass, we also require 
that each host halo have a minimum of 300 DM and 80 stellar parti¬ 
cles, which in Illustris-1 corresponds to masses of Mdm = 2 x 10^ 
M© and M* = 10^ M©. In practice, the BH mass constraint re¬ 
moves most of the lowest-mass dwarf galaxies from our sample, 
and those remaining are largely satellites that undergo minor merg¬ 
ers with larger galaxies. We find that less than 0.1% of mergers 
have a total halo mass <10^® M©, and ~ 0.1% have a total stellar 
mass < 10^ M©. 

The BH merger timescales for more massive BHs can also be 
underestimated in the simulation to a lesser degree, owing to the 
repositioning scheme and the resolution limit, which impose the 
condition that BH binary inspiral is always rapid on sub-resolution 
scales. This may be a reasonable assumption in the case of gas- 
rich major mergers, where circumnuclear gas disks and highly 
asymmetric stellar potentials can greatly reduce the BH inspiral 
times cales (e.g.,|Escala et al.|2005[[Berczik et al.|2006[|Dotti et ah] 
|2007[ [Mayer et al.|2007| l, but BH binary lifetimes in gas-poor or 

minor mergers are uncertain and may be substantially longer. 

To account for this uncertainty, we impose a delay for all BH 
mergers between the merger time in the simulation and the merger 
time assumed in our analysis. The time delay is assumed to scale 
with the BH mass ratio and the host galaxy gas fraction. Specif¬ 


ically, we add a time delay of 0.1 Gyr ! qio each merger, where 
q is the BH mass ratio, motivated by the 1/Mbh scaling of the 
[Chandrasekhar I (T94^ dynamical friction timescale for a given pri¬ 
mary galaxy mass. Because gas-driven BH inspiral should be effi¬ 
cient only in the regime where the gas disk is massive compared to 
the BHs, we add an additional 0.5 Gyr delay if the merged host 
galaxy has a cold (star-forming) gas fraction of /gas,sf < 0.1, 
where /gas,sf = Mgas,sf/(M* + Mgas,sf). Our results are not sen¬ 
sitive to the exact choice of these parameters; in practice, the effect 
of the delay imposed on gas-poor mergers is subdominant to the 
^-dependent delay. With this prescription, 26% of BH “mergers” in 
Illustris are still unmerged at z = 0, but only 2.5% of binaries with 
log q > —1.5 have not coalesced by = 0. This delayed merger 
time is used as the actual merger time in all subsequent analysis, 
and only mergers occurring at z > 0 are considered. 

In addition to the simulation snapshots, data for each BH 
merger in Illustris, as well as gas accretion data during active BH 
phases, is written at much higher time resolution|^We note that a 
small fraction of these data files for Illustris-1 in the redshift range 
0.15-0.38 were corrupted. This has only a modest effect on the 
statistics of our results, well within the uncertainties from model 
assumptions and resolution convergence. Except where otherwise 
noted, the results below refer to the Illustris-1 simulation. We dis¬ 
cuss the convergence of the Illustris-1, 2, and 3 results in Appendix 
1 ^ 

We also note that in some snapshots, a small number of black 
holes are not associated with any subhalo. Detailed analysis indi¬ 
cates that this can occur if there is a mismatch between the time 
when a subhalo is determined to have merged with a larger halo 
or EOE group and the time when its central BH is repositioned on 
the new potential minimum. Such phases are short-lived and affect 
< 0.7% of BHs at any point in the simulation. We therefore ignore 
these BHs in our analysis, with negligible impact on our results. 


2.2 BH spin models and kick velocities 

The recoil velocity of a merged BH depends on the mass ratio and 
spin vectors of the progenitor BHs. We obtain the BH mass ratio 
just prior to merger directly from the simulation data. Eor the spins, 
we must assume a physically motivated distribution. As discussed 
above, BH spins will be preferentially aligned with the orbital an¬ 
gular momentum and with each other prior to merger if their in¬ 
spiral is driven by torques from a circumbinary gas disk. If the BH 
inspiral is instead driven predominantly by stellar interactions, or 
if the accretion flow is chaotic, the spin orientations may be closer 
to random. Coherent gas accretion will also spin up the BHs on 
average, while repeated mergers at random orientations cause BH 
spin magnitudes to tend toward an approximately thermal distribu¬ 
tion. The efficiency of each of these processes in various merger 
environments is not known, but by considering a range of models 
for BH spins, we can explore their relative effects on the resulting 
recoil kicks and on the observability of offset AGN. 

The most optimistic assumption for BH spins - yielding the 
largest kicks - is that the BHs are always rapidly spinning and 
that their spins are randomly oriented prior to merger. We refer 
to this (using a dimensionless spin parameter a = 0.9) as the 


^ These data will be made publicly available in the coming months at 
the permanent site www. illustris-pro ject. org/data, where the 
other simulation data are already available jNelson et al.[2015| . 
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random-high model. We also consider a spin distribution result¬ 
ing from repeated randomly-oriented mergers, which could arise 
from a succession of dry mergers. [Lousto et ^ ( |2012| ) fit this spin 
magnitude distribution with a beta function that peaks at a ~ 0.7 
and has a large tail toward low spins (Figure [^. This is denoted as 
the random-dry model. 

At the other extreme, we may assume that BH spins always 
undergo some amount of alignment prior to merger. We use the 
hot and cold spin evolution models as described in |Dotti et al.| 
( |2010| ) and [Lousto et al.| ( [2012| ), which result in partial alignment 
(within ~30°) or near alignment ( < 10°) of BH spins, respectively. 
We adopt the analytic approximations to the spin angle and magni¬ 
tude distributions given in [Lousto et ar] ( |2012| ), which are based on 
the simulations of [Dotti et al.[ ( [2010[ ). Finally, we consider the case 
in which spin alignment is always very efficient, such that spins are 
aligned to within 5° of the binary orbital angular momentum. The 
misalignment angle of each BH is chosen at random in cos(^) over 
the interval (0°,5°), and the spin magnitude is assumed to be high 
(a = 0.9). We refer to this as the 5deg model. 

Of course, it is likely that reality lies somewhere between 
these extremes, with efficient spin alignment occurring in some BH 
mergers and not in others. Recent studies of BH spin evolution us¬ 
ing semi-analytic models have considered a mixture of coherent 
and chaotic accretion ( [Barausse|2012[[Sesana et al.|2014[ ). [Sesana[ 
[et al.[ ( [20T4] ) find that empirical constraints on BH spins are better 
fit by this type of intermediate scenario than if accretion is assumed 
to be always coherent or always chaotic. As binary BH spin align¬ 
ment may be driven by circumbinary gas disks, we construct hybrid 
spin models by utilizing information about the cold gas content of 
the progenitor host galaxies in Illustris. In general, circumnuclear 
disks with Mdisk ^ Mbh can significantly influence the inspiral- 
ing BHs. We do not attempt to estimate the fraction of gas in the 
progenitor galaxies that will end up in a circumnuclear disk, but 
rather assume that gas-rich galaxy mergers are likely to trigger the 
formation of massive circumnuclear disks. Thus, we define galax¬ 
ies as “gas-rich” if they have a cold (star-forming) gas fraction of 
> 10%, where again the gas fraction /gas,sf is defined relative to 
the total stellar mass of the progenitor subhalos. With this defini¬ 
tion, the large majority of galaxy mergers are classified as gas-rich 
- more than 85% of those at z < 1, and nearly all of those at 
higher redshift. We note that the gas-rich fraction could be over¬ 
represented in Illustris, as feedback processes do not sufficiently 
suppress star formation at the high- and low-mass end of the halo 
mass function (e.g., [Vogelsberger et al.|2014b|[Genel et al.|2014[ 
[Snyder et al.[2015[ ). However, the stellar mass function is otherwise 
in reasonable agreement with observations, as are the molecular gas 
fractions and the star formation rate density over cosmic time < |Vo-[ 
[gelsberger et al.|2014b[[Genel et al.|20T4] |. Moreover, we find that 
a lower gas-rich merger fraction would only increase the observ¬ 
ability of recoiling AGN. Alternate definitions of gas-rich mergers 
are discussed in detail in Section [377| but our results are not very 
sensitive to the exact definition used. 

We define one hybrid model in which BHs in gas-rich merg¬ 
ers have nearly-aligned spins (drawn from the cold distribu¬ 
tion), and BHs in gas-poor mergers have spins drawn from the 
random-high model. We denote this spin model as f gas-a. We 
also consider a more conservative model in which the same critical 
/gas,sf determines whether spins are drawn from the random-dry 
model (for gas-poor mergers) or the 5deg model (for gas-rich 
mergers); this hybrid model is denoted f gas-b. 

Throughout this paper, we refer to these three classes of spin 
models as “random”, “hybrid” and “aligned”. Tablesummarizes 


the definitions of each spin model. We finally note that the inclusion 
of relativistic precession effects on BH spins is beyond the scope 
of this work, but that recently developed efficient techniques for 
calculating spin precession create a promising avenue for future 
studies ( [Kesden et al.|2015|[Gerosa et al.|2015| l. 

Once the pre-merger BH spins have been obtained from the 
distribution for a given spin model, the recoil kick velocity can be 
calculated using a fitting formula based on results from numeri¬ 
cal relativity simulations. We use the folllowing formula given by 
[Lousto et al.[ ( |2012| ): 


Vrecoii = 'yme_L,i +u±(cos^e^,i +sin^e±,2) +^11^11, 


-4?7(1 + Br]), 


V± 


( 1 +^ 


(a2|| - gain), 


( 1 ) 

( 2 ) 

(3) 


«ii = (lyg + ^a5|| +145^ + X (4) 

I a 2 ± - gai±| cos(0A - 0i), 

where r] = ^/(l + ^)^ is the symmetric mass ratio, _L and || refer to 
vector components perpendicular and parallel to the orbital angular 
momentum, respectively, and and 0^,2 are orthogonal unit 
vectors in the orbital plane. The vector S = 2(a2 + g^ai) / (1 + g)^ , 
and 0A is the angle between the in-plane component A a of the 
vector A = M^(a 2 — gai)/(l + q) and the infall direction at 
merger. The phase angle depends on the initial conditions of the 
binary and is assumed to be random. The best-fit values of A = 
1.2 X 10^ km s“^, B = —0.93, H = 6.9 x 10^ km s“^, and 
^ = 145° are taken from [Gonzalez et ar] ( |2007| ) and [Lousto &[ 
[Zlochower[ p008[), and the coefficients lAg = 3677.76, Va = 
2481.21, Vb = 1792.45, and Vc = 1506.52 (all in km s“^) are 
defined in [Lousto et al.[ ( [201^ . 

The kick distributions resulting from each spin model are 
shown in Figure [^assuming a log-uniform mass ratio distribution 
in the range — 2 < log(^) < 0. It is clear that the high-velocity tail 
is greatly suppressed when spin alignment is efficient. 


2.3 Host galaxy model 

To model the merged host galaxy, we first match the BHs in each 
merging pair to their progenitor subhalos using data from the clos¬ 
est prior simulation snapshot. While the simulations provide de¬ 
tailed information about the density profiles of the progenitor host 
galaxies of merging BHs, the structure of the merger remnant at 
the time of the BH merger (which generally occurs between snap¬ 
shots) is not known. Major mergers in particular should produce 
remnant galaxies with drastically different structure than their pro¬ 
genitors, and these systems are most likely to produce observable 
recoils. Thus, we opt to use the mass and stellar half-mass radii of 
the progenitors to construct an analytic model for the gravitational 
potential of the merger remnant. 

In many cases, the BHs already inhabit a common subhalo in 
the snapshot prior to their merger, and we use the (DM, stellar, and 
gas) mass of this common subhalo as the total mass of the merged 
host galaxy. If the BHs inhabit unique subhalos, the sums of the 
DM, stellar, and gas masses of each progenitor are used for the 
total mass in the merger remnant. We note that [Rodriguez-Gomez[ 
[et al.[ ( [2015[ ) have found that the mass ratio of two merging sub¬ 
halos can vary dramatically between the time of first infall and 
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Figure 1. Normalized distributions of spin magnitude, misalignment angle, 
and recoil kick velocity. Top panel: spin magnitude distributions for the 
hot (red), cold (blue), and “dry merger” (green) models, as defined in the 
text; for comparison, the “high” spin value of a = 0.9 (black) is also shown. 
Middle panel: Distributions of the angle of misalignment between the BH 
spin vectors prior to merger, 0i^2, for the models used. The black curve 
is a distribution uniform in cos(0i,2); for clarity, we map 6 —180° — 6 
for 01,2 > 90°. The red and blue curves again denote the hot and cold 
distributions, respectively. For the orange curve, the misalignment between 
each BH spin and the orbital angular momentum (A0) is constrained to be 
< 5°; the distribution is uniform in cos(0) over this range. Bottom panel: 
recoil kick distributions resulting from Equation ^ for the spin models 
used, assuming the mass ratio q has a log-uniform distribution in the range 
—2 < log(q) < 0. The random-high model uses a = 0.9 and random 
orientations, the random-dry model uses the “dry” spin magnitude dis¬ 
tribution and random orientations, the hot and cold models correspond to 
partially- and nearly-aligned spins, and the 5deg model assumes a = 0.9 
and A6 < 5°. The high-velocity tail of kicks is strongly suppressed when 
spin alignment is efficient. 


the time of halo merger, owing to the difficulty in assigning par¬ 
ticles between two nearby halos. However, the combined mass of 
the merging halos is generally robust, and this is the quantity used 
in our calculations. This provides further motivation for our choice 
to construct analytic merger-remnant potential models rather than 
extracting density profiles directly from the simulation. Again, we 
consider only BH mergers for which each progenitor has at least 
300 DM and 80 stellar particles. 

The DM component is modeled as a |Hernquist| ( |1990| ) po- 



Gas-rich 

Gas 

-poor 


Name 

|ai,2| 

COS(0l,2) 

|ai,2| 

COS(0l,2) 

/gas,sf 

Random 

random-high 

0.9 

(-1,1) 

0.9 

(-1,1) 

- 

random-dry 

dry 

(-1,1) 

dry 

(-1,1) 

- 



Hybrid 




fgas-a 

cold 

cold 

0.9 

(-1,1) 

0.1 

fgas-b 

0.9 

(0.996,1) 

dry 

(-1,1) 

0.1 



Aligned 




hot 

hot 

hot 

hot 

hot 

- 

cold 

cold 

cold 

cold 

cold 

- 

5deg 

0.9 

(0.996,1) 

0.9 

(0.996,1) 

- 


Table 1. BH spin models used to determine the recoil kick distributions. 
We consider three classes of spin models: random, hybrid, and aligned. 
Columns: (1) Name of model, (2) spin magnitude distribution for gas-rich 
mergers, (3) distribution of spin misalignment angle (relative to the orbital 
angular momentum) for gas-rich mergers, (4) spin magnitude distribution 
for gas-poor mergers, (5) spin misalignment angle distribution for gas-poor 
mergers, (6) critical fraction of star-forming gas used to define gas-rich ver¬ 
sus gas-poor mergers (/gas,sf = ^gas,sf/(^gas,sf + M^)). The spin 
angle and magnitude distributions denoted as “dry”, “cold”, and “hot” are 
taken from |Lousto et al.|f2012) and shown in Figure[T] Otherwise, the spin 
angle distributions are randomized in cos(0) over the interval shown. 


tentia l matched to an equivalent NFW profile ( [Navarro et al.| 
\99l\ with concentration para meter Cyir — 10.5(Mvir/10^^/i 


Mo)“ (1 -Tz)~ (following Bullock et al.|2001 Maccio et al. 
|2007| ). The stellar component is modeled as a spherical bulge with 
a softened isothermal density profile, p* = cr^/(27rG(r^ + 
with velocity dispersion a‘i — GM*/i^buige and rgoft = 
where the influence radius rinfi of the BH is The pro¬ 

file is truncated at an inner radius i?ej within which matter remains 


bound to the recoiling BH. This truncated, softened profile accounts 
for the removal of the cusp of stars that is carried along with the BH 
and prevents unphysically large central densities. 

We truncate the stellar bulge at an outer radius i^buige as well; 
this radius is scaled to the stellar half-mass radii of the progenitor 
galaxies. Specifically, we take the maximum radius of the stellar 
bulge in the merger remnant to be a simple average of the progen¬ 
itor half-mass radii; this choice is based on the assumption that the 
inner density profile is steep during a major merger, and is also 
meant to account for the fact that the stellar half-light radii in II- 
lustris are up to a factor of two larger than in observed galaxies at 
z = 0, at least at the low mass end ( [Snyder et al.||2015j ). To ad¬ 
ditionally ensure that the stellar bulges do not have artificially low 
densities in low-mass galaxies, we compare the calculated stellar 
velocity dispersion to that inferred from the Mbh - relation of 
[McConnell & Ma[ ( [2013[ ). If the calculated cr* derived from M* 
and i^buige falls more than 3a below the empirical relation (assum¬ 
ing 0.4 dex scatter in log Mbh), we set cr* to lie on the relation 
and calculate i^buige from this instead of using the value inferred 
from the simulation. In practice, this criterion is imposed on most 
low-mass bulges (M* <1-2 xl0^° M©), affecting 14% of all 
merger remnants. Finally, an absolute maximum bulge radius of 15 
kpc is also imposed; this affects about 1% of all merger hosts. The 
assumption of compact stellar bulges likely overestimates the inner 
stellar density in some cases, particularly in dry mergers between 
elliptical galaxies, when the BH inspiral timescale is very long, or 
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if a disk is able to re-form after the merger. However, we consider 
this to be a conservative choice, as the steeper inner profile ensures 
that we do not overestimate the fraction of recoiling BHs that can 
escape the bulge. 

We also include a gas disk in our galaxy models; hydrody¬ 
namic simulations of recoiling BHs in isolated galaxy mergers indi¬ 
cate that the dense, circumnuclear gas disks that form during major 
gas-rich mergers can strongly suppress recoil trajectories ( |Blecha| 
|et al.||20TH ISijacki et al.||201l| |. We assume that the cold, star¬ 
forming gas in the progenitor halos has condensed into a circum¬ 
nuclear gas disk by the time of the merger, and this gas is mod¬ 
eled with a 1/r Mestel surface density profile. This is a reasonable 
approximation to the surface density profiles obtained from zoom- 
in simulations of gas-rich mergers by [Hopkins & Quataert| ( |2010] ), 
spanning radii of ~ 0.1 pc - 1 kpc. Motivated by these results, we 
assume the inner disk surface density scales with the star-forming 
gas fraction as log Eo.i = 2 log (/gas,sf/0.1) + 12 , where Eo.i is 
the surface density at 0.1 pc in M 0 /kpc^. The disk component is 
ignored in galaxies with log fgas,sf < —1.5. Typical disk masses 
are log (Mdisk/M©) ~9.5 — 10 . 5 at 2 :<l and log (Mdisk/M©) 
~10 — llatz>l. The disk is truncated at an inner radius i^ej 
within which gas is bound to the recoiling BH, and at an outer ra¬ 
dius i^max = G^Mdisk/^o» where the fiducial circular velocity vq 
of the Mestel disk is set by the above scaling of Eo.i. An absolute 
maximum disk radius of 15 kpc is also imposed; this affects about 
6 % of merger hosts. 

2.4 Recoil trajectories 

Once the merged galaxy model is assigned, the recoiling BH tra¬ 
jectory is integrated in this potential for each recoil event with 
^kAesc> 0.1. All spatial and velocity offsets are calculated relative 
to the center of the static host potential. We include stellar dynam¬ 
ical friction via the Chandrasekhar formula ( [Chandrasekhar j 1 943] ) , 
taking Mbh + Mdisk,ej as the total mass of the recoiling object. 
Gas dynamical friction is neglected; because the disk is assumed 
to be thin (/i/r < 0.1 for gas-rich systems), only a few percent of 
recoiling BHs will be kicked directly into the disk if their orienta¬ 
tions are random. If BH spin alignment is efficient and the angu¬ 
lar momenta of the accretion disk and the large-scale disk are also 
aligned, there may be a preference for kicks into the disk plane, but 
in this case the kicks will be small and difficult to observe. For large 
kicks, the v~‘^ scaling of dynamical friction ensures that it will be 
relatively inefficient as the BH leaves the center of the galaxy. The 
main uncertainty introduced by neglecting gas dynamical friction is 
its larger possible effect on recoil events with moderate Uk/Uesc that 
return to the center on short timescales, in which case the damping 
of their motion may be underestimated in our model. However, we 
demonstrate below that these objects do not dominate the popula¬ 
tion of observable recoils. 

The integration is stopped when one of the following occurs: 
a) the velocity and galactocentric distance of the BH are less than 
critical values v < 0.2cr* and R < b) the AGN luminosity 
falls below the observable limit and Uk>Uesc, or c) the integration 
reaches a Hubble time. 

Because the GW recoil kicks are implemented in post¬ 
processing, we must explicitly remove from our analysis merg¬ 
ers involving a BH that either a) was previously ejected from its 
host galaxy or b) experienced a previous recoil kick and has not 
yet returned to the galactic center. For the fiducial galaxy potential 
model, this removes 5% of BH mergers that meet all other criteria, 
if the random-high spin model is assumed. 2 % of mergers are 


excluded in the random-dry model, and < 1 % are excluded in 
all other spin models, where large kicks are less common. 

Our final sample therefore includes BH mergers that meet the 
following criteria: (i) each BH has Mbh > 10® Mq, (ii) each BH 
is hosted in a progenitor with > 300 DM (80 stellar) particles, (iii) 
the “delayed” BH merger time (defined above) occurs at z > 0, and 
(iv) the merger occurs at a sufficiently high redshift to appear in the 
past light cone of an observer at = 0 situated in the center of the 
simulation box. (Because we are concerned only with the statis¬ 
tics of observable recoils over the entire sky, this is relevant only 
at extremely low redshifts; position in the box is a negligible con¬ 
sideration at cosmological redshifts). For Illustris-1, this yields a 
sample of 8993 BH mergers. Finally, we ensure that (v) neither BH 
is still displaced from a previous recoil event at the time of merger, 
which primarily affects the random spin models. The final sample 
sizes for the random-high and random-dry spin models are 
8576 and 8800, respectively. 

2.5 AGN luminosities 

[Blecha & Loeb[ ( [2008[ ); [Blecha et al.[ ( [201 1[ ) have developed an 
analytic model for the mass and extent of the accretion disk car¬ 
ried along with a recoiling BH that is actively accreting at the 
time of the kick. The model assumes a thin “a-disk” that becomes 
self-gravitating beyond the radius where the Toomre Q parameter 
( [Toomre[[1964[ ) equals unity. The radius of the disk bound to the 
ejected BH is that at which the orbital speed Uorb ~ ^k. Typical 
ejected disk masses are a few percent of the BH mass for observ¬ 
able offset AGN (for further details see Figure]^. 

The BH accretion rate at the time of the kick is taken di¬ 
rectly from the simulation. The ejected, isolated accretion disk will 
then diffuse outward over time, causing the accretion rate to mono- 
tonically decline oc (assuming Thomson scattering opac¬ 

ity, E annizzo et al.[[1990| [Pringle[[1991|l. Following [Blecha et al.[ 
( [201 l[ l, we use this to calculate the accretion rate at each integra¬ 
tion timestep of the recoil trajectory. The accretion rate is converted 
into a bolometric AGN luminosity Lboi = eMc^, where the radia¬ 
tive efficienty e is assumed to be 0.1 at high Eddington ratios. At 
low Eddington ratios, the accretion fiow is assumed to become ra- 
diatively inefficient. Eollowing [Narayan & McClintock[ ( |2008] ), e is 
multiplied by an accretion rate dependent factor M/(/riafMEdd) 
when the Eddington ratio falls below a critical value /riaf • We adopt 
/riaf = 0.05, which is the critical Eddington ratio for the transi¬ 
tion to “radio mode” feedback in Illustris. Also, because very low 
Eddington ratios are unrealistic in the context of our thin a-disk 
accretion model, a minimum Eddington ratio must be imposed to 
prevent unfeasibly long offset AGN lifetimes. We conservatively 
adopt a minimum observable Eddington ratio of 10“^. While lower 
Eddington ratios may be detectable for some massive BHs, this 
choice is motivated by the distribution of Eddington ratios in ob¬ 
served quasars ( [Shen & Kelly[2012] ). Section [T7] examines how our 
results are affected if a lower minimum value of f^dd is assumed. 

Our procedure allows the observable offset AGN lifetime to 
be calculated for arbitrary telescope resolution and sensitivity. Eor 
each recoil event, a random viewing angle is assigned, and pro¬ 
jected separations (Ai^proj) of twice the angular resolution are de¬ 
fined as resolvable spatial offsets. Eor most of our analysis, the min¬ 
imum resolvable LOS velocity offset is assumed to be Aulos > 
600 km s“^. However, in some spectra, especially those with very 
broad lines (EWHM ^ 1000 km s“^), offsets this small may not 
be measurable. Thus, we also consider a more conservative mini¬ 
mum Aulos = 1000 km s“^. 
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The minimum observable Lboi for each event is determined by 
first defining a desired flux sensitivity and band, and then applying 
an inverse bolometric correction to the corresponding luminosity at 
the redshift of the recoil event. We use the luminosity-dependent 
optical and hard X-ray bolometric corrections from [Hopkins et al.] 
( |2007| ). For the near-IR, we assume a constant bolometric correc¬ 
tion of 11.0, derived from the mean SED of [Richards et al.| ( [2006| ) 
at 2 /xm. AGN bolometric corrections are not well constrained, par¬ 
ticularly at optical wavelengths where obscuration and host galaxy 
contamination have large effects. Even at X-ray wavelengths, the 
fraction of Compton-thick sources is uncertain, and bolometric cor¬ 
rections have been found to vary with Eddington ratio ( [Vasudev^ 
[& Eabian|2007[[Lusso et al.|2012[ l. We neglect obscuration effects 
here, as these should be less significant for ojfset AGN (though the 
early post-recoil stage of velocity-offset AGN may be a notable ex¬ 
ception, as discussed in Section [T^ . 

Despite the inherent uncertainties, our results depend only 
weakly on the choice of bolometric correction. This is primarily 
because at low to moderate redshift, the minimum f^dd = 10“^ 
criterion is typically a more stringent constraint on the AGN life¬ 
time than the absolute flux limit derived from the (reverse) bolo¬ 
metric correction. The consequences of this are discussed in detail 
in Section [331 


3 RESULTS 

3.1 Merging BH population 

The differential all-sky merger rate of BHs in Illustris at z = 0, 
d^A^/dlog(;2;) dt, is shown in Eigure|^ Two BHs are allowed to 
merge in the simulation when their separation is less than a smooth¬ 
ing length. In cosmological simulations, this spatial scale (~ 0.1-1 
kpc in Illustris) is far larger than the regime in which BH inspiral 
is dominated by GW emission. Thus, as described above, we add 
a delay timescale tdeiay = O.lGyr/^ + tgas to the time of each 
BH merger, where tgas = 0.5 Gyr if the fraction of mass in cold 
(star-forming) gas /gas,sf <0.1 and zero otherwise. More detailed 
modeling of binary BH inspiral on sub-resolution scales is beyond 
the scope of this work, and will be examined in an upcoming study 
(Kelley et al., in preparation). 

The green lines in Eigure|^show the merger rate in different 
mass bins, for the remnant BH mass Mtot = Mi + M2. Over a 
large redshift range 0.5 < z < 3, the merger rate is dominated by 
BHs in the mass range 10^-10® M©. Lower-mass BHs dominate at 
higher redshifts, and higher-mass BHs dominate at lower redshifts, 
consistent with expectations from hierarchical growth. Low-mass 
BH mergers (^ 10® M©) are of interest as GW sources for a future 
space-based interferometer such as eLISA, but these are excluded 
from our analysis owing to resolution limits, as discussed above. 

The cumulative merger rate of this BH sample at z = 0, inte¬ 
grated over cosmic time, is 0.42 yr“^, and the merger rate for BHs 
with total mass 10^ - 10® M© is 0.27 yr“^. This is in reasonable 
agreement with results from semi-analytic modeling (e.g., [Sesana[ 
[et al.[[2004| [Barausse[[2012[ l, especially considering the very dif¬ 
ferent BH seeding and growth prescriptions in the semi-analytic 
models versus Illustris. The BH merger rate also peaks at simi¬ 
lar redshifts {z > 2) in both Illustris and the semi-analytic models, 
though after the merger delay timescale is implemented, the peak 
shifts to slightly lower redshift {z ~ 1.5) in our sample. However, 
there are more high-redshift mergers in the semi-analytic models, 
owing to the BH seed formation prescriptions, which form seeds at 



Figure 2. The BH merger rate in Illustris at 2 ; = 0, per unit redshift (in 
logarithmic bins). The thick gray curve denotes the “raw” merger rate for 
all BH mergers that meet the selection criteria defined above, except that 
the actual BH merger redshift from the simulation is assumed instead of the 
“delayed” redshift. The thick black curve shows the merger rate with a delay 
timescale added to each merger that depends on BH mass ratio and host gas 
fraction, as described in the text. The thin green lines denote the merger rate 
for different remnant BH mass bins, as indicated in the plot legend. The raw 
merger rate peaks at 2 ; ~ 2, in good agreement with semi-analytic models, 
and the delayed merger rate peaks at 2 ; ~ 1.5. The cumulative merger rate 
over cosmic time also agrees well with semi-analytic models, but more low- 
redshift mergers occur in Illustris. 


higher redshifts than in Illustris. The first BHs in Illustris form at 
2: ~ 10 - 11, and the first mergers occur at z < 7. The merger rate 
at low redshift {z < 1) is correspondingly higher in Illustris than 
in the models of [Sesana et al.[ ( |T004[ ) (for BH binaries with tota l 
mass >10® M©). In the “heavy seed” model of 


Barausse 


( 20121 , 


the cumulative merger rate (for Mtot >10® M©) is about a fac¬ 
tor of three higher than in Illustris, such that the merger rate per 
unit redshift is slightly higher than the Illustris rate at all redshifts. 
Eurther comparison with comparable cuts on host halo masses im¬ 
proves the agreement with both the [Baraussej \2012\ and [Sesana[ 
[et al.[ ( |2004) models, though some differences in the mass and red¬ 
shift distribution of the merger rates remain, as does the overall 
trend toward more high-redshift mergers in the semi-analytic mod¬ 
els (E. Barausse and A. Sesana, private communication). Empiri¬ 
cal constraints cannot yet distinguish between these differences in 
the details of BH formation, growth, and merger histories. These 
distinctions between cosmological simulations and semi-analytic 
models are important to bear in mind, and they contribute to the 
differences in offset AGN rates predicted here and in [Volonteri &[ 
[Madau[ ( |2008[ ) (Section [43] ). However, the level of agreement of the 
overall merger rates in Illustris versus semi-analytic models is en¬ 
couraging. We finally note that the galaxy merger rates in Illustris 
are in reasonable agreement with observations ( [Rodriguez-Gomez[ 
[et al.|2015] ), though they disagree qualitatively with rates derived 
from semi-analytic models based on cosmological N-body simula¬ 
tions <Guo & White|2008> . 
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Figure 3. The merging BH mass ratio distribution in Illustris, separated into 
total BH mass bins. As indicated in the plot legend, the blue curve shows 
the distribution for all BH masses (above the minimum 2 x 10® Mq, and 
the cyan, green, and orange curves correspond to remnant BH masses of 
log (Mbh/Mq) > 7, 8, and 9, respectively. For most BH masses, the q 
distribution is nearly flat in log space over the range -2 < log < 0. The 
dotted lines show the distributions for only BH mergers occurring at 2 : < 1. 
The thin magenta line denotes mergers that yield a recoil velocity > 0.1 
Vesc', no mergers with log g < — 1.3 produce signiflcant recoils. 

The most massive BH mergers (Mtot ^10^ M©) will create 
luminous GW signals that could be detected with PTAs; we find a 
cumulative merger rate of 0.006 yr“^ for these massive mergers. 
This is consistent with predictions from semi-analytic models that 
individual GW events resolvable with PTAs should be rare (e.g., 
[Sesana et al.|2009| l. However, the stochastic background from un¬ 
resolved massive BH mergers may be detectable with PTAs in the 
coming years. The predicted GW signal from BH mergers in Illus¬ 
tris will be explored in detail in Kelley et al. (in prep). 

Figure shows the distribution of merging BH mass ratios q, 
separated by total BH binary mass. For a wide range of BH masses, 
the log q distribution is nearly fiat over the range — 2 < log ^ < 0. 
Low mass ratios are favored slightly for total BH masses > 10® 
M©; these massive BHs are fewer in number, so they more often 
undergo unequal-mass mergers. The trend reverses when mergers 
with total mass Mbh = 10^ - 10® M© are included, and when all 
BH masses >10® M© are considered, the distribution clearly rises 
toward log q = 0. Recall that a lower mass limit of 10® M© has 
been imposed for individual merging BHs, such that BH binaries 
with total mass < 10^ M© cannot have log q < —1. 

The dotted lines indicate the q distribution for only mergers 
occurring at z < 1. We see that major mergers are more common 
at high redshift. This is partly a result of hierarchical growth; minor 
mergers become more frequent as larger galaxies (and BHs) form. 
It is also a product of the long dynamical friction timescales for 
minor mergers, refiected in the ^-dependent delay timescales we 
have imposed on the simulation merger time. The limited box size 
of Illustris, (106.5 Mpc)®, could further contribute to the limited 
number of major mergers between low-redshift, massive BHs. For 
mergers at z < 1, the distribution is nearly fiat or declining with 
log q for all BH masses. 

Also shown in this figure is the subset of BH mergers that pro¬ 
duce a recoil kick with v^Vesc >0.1 (in the random-dry spin 
model), which is the minimum kick velocity for which recoil tra¬ 


jectories are calculated. Only mergers with log g > — 1.5 produce 
such recoils, with most requiring log ^ > — 1. 

We note that the merging BH mass ratio distribution in Illustris 
is qualitatively different from that in the models of VM08, which 
exhibits a sharp cutoff at g > 0.3. Because the recoil kick distribu¬ 
tion depends sensitively on BH mass ratios, this is an important dis¬ 
tinction. While there are uncertainties in the BH merger timescales 
in Illustris, as discussed above, they are largest for minor merg¬ 
ers. A more detailed comparison of our results with those of VM08 
is given in Section |4.2| We additionally note that recent numeri¬ 
cal simulations of BH binaries in circumnuclear gas disks suggest 
that accretion onto the lower-mass BH may be more rapid, driving 
the final mass ratios closer to unity than the progenitor distribution 
( |Roedig et al.|20 1 2| [Farris et al.|2014| ). 


3.2 Recoiling BH dynamics 

Figure]^ shows the distribution of kick velocities scaled to the host 
escape speed, for a random spin model (random-dry), a hybrid 
model (fgas-b), and an aligned model (5 deg). For the random 
spin models, the distribution peaks at v^Vesc ~ 0.3, unlike the 
hybrid and aligned models, for which the distribution is heavily 
skewed toward low Uk/Uesc- In the 5 deg model, recoiling BHs 
almost never escape the galaxy entirely. (No recoils have v^Vesc 
> 1 in the 5deg model, but this is limited by mass resolution 
and by statistics. The largest kicks produced in the 5 deg model 
could escape the lowest-mass galaxies in our sample, and would 
easily escape dwarf galaxies unresolved in the simulation.) Escap¬ 
ing BHs are also very rare in the cold spin model, comprising 
< 0.05% of all kicks. The escape fractions in the hybrid spin mod¬ 
els are 0.4% and 0.2% for fgas-a and fgas-b, respectively, 
comparable to the hot spin model (0.6% escaping). In contrast, 
the random spin models eject a full 6% (random-dry) and 12% 
(random-high) of all merged BHs. 

The colored lines in Figure]^ show the distribution of v^Vesc 
separated by total merging BH mass. Because BH mass correlates 
with galaxy mass and Uk depends only on the BH mass ratio, large 
Vk/vesc events are rarer for high-mass BHs. Despite this, in the 
random and the hybrid spin models, even the most massive BHs 
(> 10^ M©) occasionally escape their hosts. Thus, if superkicks 
do occur, even massive hosts may be left without a central BH. 

For bound recoiling BHs (v^Vesc < 1), trajectories can be 
characterized by the maximum apocentric distance from the host 
(^max) and the time spent off-center before returning to the galac¬ 
tic nucleus (tretum). Thcsc quantities are plotted versus v^Vesc for 
selected spin models in Figure Naturally, both R max and treturn 
increase with Uk/Uesc, but the scatter is large, owing to the large 
dynamic range in halo, stellar, and gas model parameters. There is 
little to no correlation below v^Vesc^ 0.5, where recoil trajecto¬ 
ries are strongly suppressed by the dense stellar and gas compo¬ 
nents. These low-velocity recoiling BHs are generally confined to 
the central kpc of the galaxy, with return time < 10^ yr, and some 
are displaced by less than a parsec. In contrast, large kicks with 
Vk/vesc ^0.8 can displace the BH for > 1 Gyr, with i^max ^ 100 
kpc. 

Pre-merger spin alignment greatly suppresses high-velocity 
recoils. The hybrid fgas-b spin model produces a much smaller 
fraction of recoils with v^Vesc >0.6 than the random spin model, 
and the 5deg model yields no kicks above v^Vesc = 0.8. In the 
latter model, with extremely efficient spin alignment, recoiling BHs 
are nearly always confined to the central kpc of the galaxy (though 
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Figure 4. The distribution of kick velocities scaled to the host escape speed, V]^lvesc, is shown for a random spin model (random-dry, left panel), a hybrid 
model (fgas-b, middle panel), and an aligned model (5deg, right panel). The distributions are separated by total merging BH mass using the same color 
scheme as in Figurej^ and as indicated in the plot legend. The random model has a large tail of escaping BHs, the hybrid model has only a few, and the aligned 
model has none. 
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Figure 5. The characteristics of bound recoil trajectories (uk/r’esc < 1) are shown for three different spin models: random-dry {top panels), fgas-b 
{middle panels), and 5 deg {bottom panels). In each case, the left panel shows the distribution of return times fretum for the recoiling BH to settle back to the 
galactic center, as a function of Uk/^’esc• The right panel shows the maximum galactocentric distance of the BH orbit, i7max, versus Uk/r’esc- There is little 
dependence of either quantity on V]^/vesc for values < 0.5, but marginally bound recoils can wander in the halo for many Gyr. Such events are much rarer in 
the hybrid model than the random model, and almost never occur in the aligned spin model. 
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Figure 6. Top panel: the distribution of Eddington ratios (/Edd = 
^AGN/^Edd) is shown for merging BHs in Illustris, separated by red- 
shift bin according to the plot legend (thicker lines correspond to lower z). 
Bottom panel: same distribution, but separated by the mass fraction (rela¬ 
tive to stars) of cold, star-forming gas in the host (thicker lines correspond 
to lower /gas,sf )• Unlike the /Edd distribution for all BHs (cf.|Sijacki et al.| 
|2015) , the distribution for merging BHs heavily favors values near unity at 
all redshifts and for a wide range in gas fractions. 

again, this does not apply to dwarf galaxies that are below the res¬ 
olution limit of the simulation). 

3.3 Merging and recoiling AGN properties 

The AGN Eddington ratios (/Edd = Lagisi/L^ dd) in Illustris 
evolve strongly toward lower values with cosmic time ( [Sijacki et al.| 
|2015| Figure 9). At z = 4, when BHs have low masses and their 
hosts are generally gas-rich, the /Edd distribution peaks near unity. 
At z = 0, Eddington-limited AGN are very rare, and the median 
/Edd for all BHs is 3.7 X 10“^. 

The Eddington ratios for merging BHs are quite different. We 
find that the /Edd distribution is heavily skewed toward unity for 
merging BHs at all redshifts. Figure (top panel) shows the /Edd 
distribution in several redshift bins, for all merging BHs that meet 
our minimum-mass criteria. At low redshifts, a tail extends to very 
low Eddington ratios (/Edd < 10“^), but more than half of all 
merging BHs have accretion rates capped at the Eddington limit at 
the time of merger. The weak evolution of /Edd with redshift rel¬ 
ative to the overall BH population suggests that merger dynamics, 
rather than cosmic epoch, determine the fueling rate of these AGN. 



Figure 7. The mass distribution of the accretion disk bound to the recoiling 
BH, scaled to BH mass. The thick black line shows the distribution for all 
mergers, while the thin lines denote the distributions for only recoil events 
that produce observable spatially-offset (solid blue line) and velocity-offset 
(dashed red) AGN. For the offset AGN, the random-dry spin model 
and HST-COSMOS sensitivity (AB(F814W) = 27.2 mag for point sources, 
[Koekemoer et al.|2007t , resolution (0.1”), and Aulos > 000 km s“^ are 
assumed, and a minimum offset lifetime of 10^ yr is imposed. Typical disk 
masses are a few percent of the BH mass. 

The bottom panel of Figurej^shows the /Edd distribution sep¬ 
arated by the star-forming gas fraction in the host, illustrating that 
higher-/gas galaxies have higher Eddington ratios on average, but 
even mergers with /gas,sf <0.1 have an /Edd distribution skewed 
toward unity. This indicates that the cold gas in these galaxies is ef¬ 
ficiently funneled to the central regions during major mergers. The 
ability of major mergers to trigger AGN fueling is well-studied in 
simulations of isolated galaxy mergers (e. g.,|Wyithe & Loeb|2003[ 
|Di Matteo et al. 12005}[Hopkins et al.|2(j06l|2008at , and a clear cor- 

relation between merging activity and AGN has been observed in 


real systems (e.g., Sanders et al.||1988| |Urrutia et al.||2008| Koss 

et al. 

2010||Ellison et al.|2011||Satyapal et al.|2014| Treister et al. 

2012 

1. Our results demonstrate that this effect is self-consistently 


produced in the Illustris cosmological simulations, and that a strong 
correlation exists during the late stages of the galaxy merger. In 
cases where a significant delay occurs between the BH merger time 
in the simulation and the actual BH merger on sub-grid scales, the 
accretion rate may be lower at the time of the merger and recoil, 
but the longest delays occur for minor mergers that are less likely to 
produce observable recoils. Once the BH is ejected from the galaxy, 
however, its accretion rate declines with time, so the recoiling AGN 
have L ~ L^dd only briefly. 

Note that even though the recoiling BH accretion rate in our 
models declines monotonically after the recoil kick, the total AGN 
lifetime is generally longer than if a constant accretion rate (and 
corresponding AGN lifetime Mdisk/MeH) is assumed. For the 
minimum observable value of /Edd = 10“^ assumed in our mod¬ 
els, the AGN lifetime may be nearly 10 times longer than a constant 
Mdisk/MsH lifetime. 

As described in Section [O] recoiling BHs should carry along 
material that orbits the BH with Uorb ^^k- The amount of mass 
bound to the recoiling BH, along with the accretion rate at the 
time of the kick, are used to determine its AGN lifetime. We cal¬ 
culate the mass and radius of this bound material assuming it is 
composed of an alpha-disk that becomes self-gravitating beyond 
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Figure 8. The distributions of observable offset AGN lifetime versus 
are shown for resolvable projected spatial offsets (left panels) and 
for LOS velocity offsets (right panels), weighted by the total number of 
recoil events. HST-COSMOS sensitivity and resolution, and a minimum 
^'i^LOS > 600 km s“^, are assumed. The top panels correspond to the 
random-dry spin model, which has 8800 total recoil events, 492 of 
which yield an observable spatially-offset AGN, and 413 of which produce 
an observable velocity-offset AGN. The bottom panels correspond to the 
fgas-b model, which has 8989 total events, 54 (30) of which produce ob¬ 
servable spatial (velocity) offsets. In all cases, the longest lifetimes occur 
for kick speeds near the host escape speed. 


the radius where the Toomre Q = 1. The resulting ejected disk 
masses for recoils with v^Vesc > 0.1 are shown in Figurej?] for the 
random-dry spin model. The overall distribution of disk masses 
peaks at ~ 0.1 Mbh and has a tail extending to < 10 “^Mbh. 
Among observable offset AGN, very few have Mdisk ~ Mbh, be¬ 
cause these correspond to low-velocity recoils; this is true regard¬ 
less of the degree of spin alignment. The disk masses for observable 
velocity-offset AGN are slightly smaller than for spatially-offset 
AGN, because the minimum velocity offset requirement means that 
the former have higher kick speeds on average. In the random spin 
models, most offset AGN have Mdisk/MBH = 1 - 10%. In the hy¬ 
brid and aligned models (with lower kick velocities), typical disk 
masses are slightly higher, ~ 3 - 30% Mbh. 

Figure shows the relationship between v^Vesc and the 
spatially- or velocity-offset AGN lifetime. Here and throughout 
our analysis, a random viewing angle is assigned to each recoil 
event, and the offset lifetime is calculated as the total time for 
which the recoiling AGN is detectable and has a resolvable pro¬ 
jected spatial offset or a resolvable line-of-sight (LOS) velocity 
offset. The random-dry and fgas-b spin models are shown, 
assuming HST-COSMOS sensitivity (AB(F814W) = 27.2 mag 


for point sources, [Koekemoer et al.||2007| l, resolution (0.1”), and 
Aulos > 600 km s“\ 

The spatially-offset AGN lifetimes are longest for recoils at 
or somewhat below the escape speed (uk/uesc ~ 0.6 - 1). This is 
because marginally-bound BHs experience little deceleration and 
spend most of their time at large apocenters, well separated from 
the host centroid. Escaping BHs also produce large offsets, of 
course, but this is countered by the fact that the accretion disk mass 
carried with the BH decreases oc such that AGN lifetimes are 
vanishingly small for Uk ^ Uesc- Such extreme kicks are also rare. 
At high redshift, escaping BHs do begin to dominate the popula¬ 
tion of observable offset AGN, owing to both angular resolution 
limits and the suppression of low-velocity recoils in gas-rich merg¬ 
ers. The fact that offset AGN lifetimes are maximized for Uk ^ Vesc 
has important consequences for the observability of recoils, fore¬ 
most that spatial resolution should not be a strongly limiting factor 
in systematic searches for recoiling AGN. 

Similarly, velocity-offset AGN have the longest lifetimes for 
recoils near the escape speed (Figure right panels). On aver¬ 
age, velocity-offset AGN have shorter lifetimes than spatially offset 
AGN, because many recoil events have Uk > 600 km s“^ initially 
but are quickly decelerated in the host potential. This is especially 
true at high redshift, where gas-rich mergers tend to suppress long- 
lived recoils. 

These trends are qualitatively similar for models that include 
pre-merger BH spin alignment, except that there are fewer offset 
AGN overall, and the high v^Vesc tail of events is mostly trun¬ 
cated (Figurebottom panels). Recoil events in gas-rich mergers 
that produce short-lived velocity-offset AGN in the random-dry 
model are mostly suppressed in the hybrid model. 


3.4 Spatial and velocity offsets 

Figure shows time-weighted distributions of projected spatial 
(Ai7proj) and velocity (Aulos) offsets for all recoiling AGN in the 
simulation (i.e., not only those with resolvable offsets are shown). 
The panels compare the distributions for the random-dry, 
fgas-b, and 5deg spin models. For the random-dry model, 
the distribution is dominated by recoiling AGN with velocities 
above a few hundred km s“^ and spatial separations ^0.1” (physi¬ 
cal separations > 1 kpc are common). This follows from the result 
discussed above: recoils with Uk ^ Vesc have the longest offset 
lifetimes. The largest offsets (> 1000 km s“^ and > 10 kpc) re¬ 
sult primarily from BHs that escape their host galaxies, and the an¬ 
gular offset distribution has a tail that extends beyond 10”. This 
demonstrates that seeing-limited surveys could resolve spatially- 
offset AGN. 

The right panel of Figure demonstrates that if spins are al¬ 
ways nearly aligned, spatial offsets > 0.1” are rare, and resolvable 
velocity offsets (> 600 km s“^) never occur in the 5deg spin 
model. The confirmation of a single recoiling AGN with a resolv¬ 
able velocity offset would therefore be a strong indication of spin 
misalignment ^5°. We also see that the Ai^proj, Aulos distri¬ 
bution is no longer skewed toward the largest offsets in this case, 
owing to the lack of escaping BHs and short return times in this 
spin model. 

For the fgas-b hybrid spin model (middle panels of Figure 
1^, the distributions of Ai^proj and Aulos offsets are bimodal. 
Offsets greater than 0.5” and 600 km s“^ are produced almost 
exclusively in gas-poor mergers, for which BH spins are assumed 
to be random. Small offsets are produced mainly by mergers with 
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Figure 9. The time-weighted distributions of projected spatial offset versus LOS velocity offset are shown for all recoiling AGN. The color scale indicates 
the amount of time a BH spends at a given (Ai?proj, Aulos) while active as an AGN and offset from the host nucleus, integrated over all merger events 
over cosmic time, regardless of whether the recoil event produces a resolvable offset. The sensitivity of HST-COSMOS is assumed. The top panels show the 
distribution of angular offsets, and the bottom panels show the physical spatial offsets. The left panels correspond to the random-dry spin model, middle 
panels show the f gas-b model, and right panels show the 5 deg model. Large offsets > 1 kpc and > 1000 km s“^ are, favored in the random spin model 
but do not occur at all if BH spins are nearly aligned. 


/gas,sf > 0.1, for which BH spins are nearly aligned. This cre¬ 
ates distinct signatures in the host galaxy population relative to the 
non-hybrid spin models, as we show in Section [3^ 

In Figure the time-weighted distributions of bolometric 
AGN flux and Eddington ratio are shown versus Ai^proj and 
Aulos, for the random-dry spin model. As in Figureall re¬ 
coiling AGN are shown, not just those with resolvable offsets. We 
see that for angular separations ^1”, recoiling AGN are generally 
above the flux limit. This reflects the fact that the largest angular 
separations occur at low redshift, where the minimum AGN lumi¬ 
nosity is generally determined by the criterion f^dd > 10"^ in- 
Stead of by the flux limit. Even in the case of velocity-offset AGN, 
where there is no dependence on angular scale, offset AGN life¬ 
times at low z are usually limited by the f^dd criterion rather than 
the flux limit, at least for HST sensitivity. The bottom panels in Fig¬ 
ure [^confirm that the population of recoiling AGN is dominated 
by low Eddington ratios (< 0.1), particularly for large Ai^proj 
and Aulos offsets. (Lower-velocity recoils typically return to the 
galactic center before their luminosity can decline substantially.) 

The same distributions are shown for the 5deg spin model in 
Figure [TT] Owing to the short return times of recoils in this model, 
the flux limit is rarely reached, and the /Edd distribution is heavily 
skewed toward high values. 

One important caveat to these flndings is that our models in¬ 
clude only accretion onto recoiling BHs from the disk carried along 
with the BH at the time of the kick. As shown in IBlecha et al.l 
( |2011| ), recoiling BHs on bound orbits may encounter fresh fuel on 
subsequent passages through the galactic center, substantially ex¬ 


tending the AGN lifetime in some cases. This may create a popula¬ 
tion of offset AGN with smaller spatial offsets, such that the actual 
distribution is less biased toward large separations than our results 
suggest. We also assume the recoils always occur in a dense stellar 
cusp that forms during the galaxy merger (e.g., |Mihos & Hernquist| 
1 1994 1 [Hopkins et al.|2008b|[2009| ). If a recoil event instead occurs 
in a cored elliptical (perhaps following a dry merger, where the BH 
inspiral time may be long), the BH may have a long phase of small- 
scale oscillations about the galactic center ( jGualandris & Merrittj 
[2008[[Lena et al.|2014j ). However, any new gas supply encountered 
by the recoiling BH on small scales would only add to the offset 
lifetimes predicted by our models. 

Offset AGN that exhibit both velocity and spatial offsets, such 
as CID-42 ( [Civano et al.[2010[[2012[ ), will be especially compelling 
recoil candidates. We And that a majority of velocity-offset AGN 
(weighted by offset lifetime) will have projected spatial offsets re¬ 
solvable with high-resolution imaging (e.g., HST). This fraction in¬ 
creases with increasing spin alignment, from > 50% in the random 
spin models to ~ 90% in the cold model. Thus, high-resolution 
imaging may be a fruitful means of follow-up for any velocity- 
offset recoil candidates. Even with SDSS resolution, ~ 4 - 20% 
of velocity-offset AGN may have resolvable spatial offsets; these 
primarily occur at low redshift {z < 0.5). 

In contrast, only ~ 10 - 20% of spatially-offset AGN should 
have discernible velocity offsets Aulos > 600km s“^ (< 15% in 
the aligned spin models). With a more stringent velocity offset cri¬ 
terion (Aulos > 1000km s“^), < 10% of spatially-offset AGN 
exhibit simultaneous velocity offsets, and no spatial/velocity offset 
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Figure 10. Similar to Figure]^ the time-weighted distributions of pro¬ 
jected spatial offsets {leftpanels) and LOS velocity offsets {rightpanels) are 
shown for all recoiling AGN in the random-dry spin model (and HST- 
COSMOS sensitivity), for all timesteps during each recoil event and for all 
redshifts. Here the offsets are plotted versus the bolometric AGN luminos¬ 
ity {top panels), bolometric AGN flux {middle panels) and Eddington ratio 
^boi/^Edd {bottom panels). The Eddington ratio is constrained to have a 
maxmium of 1 and a minimum of 10“^ in our models. The feature apparent 
at /Edd = 0.05 arises from the transition to the radiatively-inefflcient ac¬ 
cretion regime. Most offset AGN have /sdd < 0.1, and for spatial offsets 
> 1”, observable offset AGN are rarely flux-limited. 


overlap occurs in the aligned spin models. However, we demon¬ 
strate below that even this small fraction of spatially-selected offset 
AGN exhibiting velocity shifts could yield a population of strong 
recoil candidates. 

One difficulty with confirming the nature of objects like CID- 
42 is that the offset AGN is superimposed on the host galaxy bulge, 
such that it cannot be easily distinguished from a member of an 
inspiraling BH pair. CID-42 has no evidence for a second AGN in 
the host nucleus, but excluding the presence of a quiescent BH is 
more difficult. The strongest offset AGN candidates will likely be 
those that have not only a resolvable spatial offset (and possibly a 
simultaneous velocity offset), but are also well-separated from the 
host galaxy’s stellar light. 

We find that a non-negligible (time-weighted) fraction of 
spatially-offset AGN resolvable with HST also have separations 
Ai^proj > ^buige, where i^buige is the radius of the stellar bulge 


Figure 11. The same distributions as in Figure are shown, but for the 
5deg spin model. The short return times of recoiling BHs in this model 
prevent the flux limit from being reached in most cases. 

component in our models. Specifically, ~ 20 - 30% of spatially 
offset AGN have Ai^proj > i^buige in the random and hybrid spin 
models, and a similar fraction of AGN with simultaneous spatial 
and velocity offsets are also found outside the stellar bulge. Separa¬ 
tions > i^buige are rare in the aligned spin models. For offset AGN 
resolvable with SDSS, ~ 40 - 50% of spatially-offset AGN may 
have Ai?proj > ^buige (in the random and hybrid spin models), 
and ~ 20 - 40% of simultaneously spatially- and velocity-offset 
AGN should be found at such separations. We caution that the stel¬ 
lar bulges in our galaxy potential models may be more compact 
than in reality if there is a long delay between the galaxy coales¬ 
cence and the BH merger and recoil, so i^buige niay be underes¬ 
timated in some cases. Nonetheless, these results indicate that a 
subsample of offset AGN displaced beyond the host stellar bulge 
could be found, and these would be ideal targets for follow-up and 
confirmation of recoiling BHs. 

3.5 Recoiling AGN observability 

3.5.1 Redshift distribution and general trends 

The predicted number of observable offset AGN per square degree 
is shown is Figureas a function of redshift, for various surveys 
and for three different spin models (random-dry, f gas-b, and 
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cold). In the random-dry spin model, 2 - 3 per deg^ spatially 
offset AGN could be detected in HST-COSMOS ( [Koekemoer et al.| 
2007| ) out to z = 3, and ~ 1 per deg^ could be found at z < 1. 
The predicted source counts decrease with survey sensitivity and 
spatial resolution; SDSS could see about 1 spatially-offset AGN 
per 40 deg^, and none should be detectable at z > 0.8. 

HST could detect velocity-offsoX AGN in slightly lower num¬ 
bers; > 1 per deg^ are predicted out to z = 3 for a minimum re¬ 
solvable offset of Aulos > 600 km s“^. In general, for angular 
resolutions better than ~ 0.5”, we find that spatially-offset AGN 
should be found in greater numbers than velocity-offset AGN. The 
reverse is true for seeing-limited surveys with spectroscopic cov¬ 
erage; in the random spin model, SDSS should find ~ 10 times 
more velocity-offset AGN than spatially-offset AGN for a mini¬ 
mum Aulos > 600 km s“^. Even if only offsets Au > 1000 
km s“^ can be resolved, velocity offsets could still be detected in 
greater numbers with SDSS. 

The predicted number of observable velocity offset AGN (for 
Aulos > 600 km s“^) differs by at most a factor of a few between 
HST and SDSS, and even less at low redshift, despite the large dif¬ 
ference in sensitivity. As discussed above, the imposed minimum 
AGN Eddington ratio of 10“^ generally corresponds to a higher 
luminosity limit than does the assumed flux limit, for low to moder¬ 
ate redshifts. Thus, our results depend only weakly on the absolute 
fiux limit, which becomes important only at higher z. 

This trend is even more prominent in the hybrid spin models, 
where there is almost no dependence on the assumed fiux limit. 
Relative to the random spin models, the hybrid models produce 
fewer superkicks, and those that do occur are limited to dry merg¬ 
ers, which are found predominantly at low redshifts (virtually none 
occur at z > 1). Thus, the AGN lifetime is almost always set by the 
minimum /Edd, and increasing the angular resolution only moder¬ 
ately increases the number of observable spatial offsets for these 
low-redshift sources (~ factor of 10-12 higher number counts for 
a factor of 20 increase in spatial resolution). Note that the SDSS 
predictions for spatially offset AGN are similar between the ran¬ 
dom and hybrid spin models; most recoil events in the random spin 
model that produce a spatial offset observable with SDSS also pro¬ 
duce an observable offset in the hybrid spin model. This refiects the 
fact that SDSS is most sensitive to luminous AGN at low to moder¬ 
ate redshift, which are more likely to occur in dry mergers that can 
produce superkicks in the hybrid model. 

In the aligned (cold) spin model, fewer observable recoils 
are predicted than in the random or hybrid models, as expected, 
and there is a steeper dependence on angular resolution. Only 1 
spatially offset AGN per 30 deg^ should be detectable with HST, 
and velocity offset AGN should occur at rate of < 1 per 100 deg^. 

We note that, despite the prevalence of major mergers among 
BHs with Mbh < 10^ Mq (Eigure |^, these low-mass BHs 
do not dominate the number counts of observable offset AGN in 
any of our models. In fact, for the hybrid spin models, massive 
BHs (Mbh > 10® M©) dominate the sample of offset AGN, 
because these BHs are preferentially hosted in gas-poor galaxies 
where superkicks can occur. Additionally, for seeing-limited sur¬ 
veys, spatially-offset AGN tend to have high BH masses regard¬ 
less of spin model. The offset AGN lifetime depends on BH mass, 
such that low-mass BHs will often exhaust their fuel supply before 
reaching a resolvable separation from the host galaxy. 

This trend occurs even though massive BHs are generally 
hosted in massive galaxies, which have high central escape speeds. 
As a result, only recoiling AGN with Uk/Uesc ^0-8 produce spatial 
offsets that are resolvable with seeing-limited surveys. One caveat 


to this is that our sample of dwarf galaxies is incomplete owing to 
mass resolution limits. A more complete sample of dwarf galax¬ 
ies could produce a population of lower-mass AGN with resolvable 
offsets. 

As mentioned above, our accretion disk model generally re¬ 
sults in longer recoiling AGN lifetimes than if a constant accre¬ 
tion rate is assumed. This is especially important for spatially-offset 
AGN, for which offset AGN lifetimes are more often limited by de¬ 
clining luminosity than by a return to the galactic center. (Velocity- 
offset AGN have shorter lifetimes on average, owing to decelera¬ 
tion by the host galaxy potential.) In most cases, our models predict 
lower rates of observable spatially-offset AGN by up to a factor of 
a few if a constant Mdisk/MBH lifetime is assumed. The difference 
is largest for the aligned spin models and seeing-limited resolution, 
where in the constant-M model, few BHs achieve a resolvable sep¬ 
aration during their active lifetime. In contrast, velocity offset AGN 
lifetimes are often limited by deceleration in the host potential, 
rather than by the luminosity limit, and the predicted numbers of 
observable velocity offset AGN can be a factor of a few lower or 
higher when a constant Mbh is assumed. 


3.5.2 Imaging surveys 

Eigure shows the cumulative number of spatially-offset AGN 
per square degree at z < 3 observable in selected current and future 
surveys, for all spin models considered in this study. We first exam¬ 
ine the random spin models: random-high and random-dry 
(left panels). Up to 7 per deg^ spatially offset AGN may be de¬ 
tectable with JWST. In the hard X-ray band, up to ~ I per deg^ off¬ 
set AGN could be detected at the sensitivity of the Chandra Deep 
Eield North (CDE-N, [Alexander et al.|206^ survey, and ~ 0.3 per 
deg^ at the limit of the Chandra-COSMOS Legacy (C-COSMOS, 
jCivano et al.|2013| l survey. Seeing-limited surveys (assuming 1” ef¬ 
fective resolution) may resolve up to 1 spatially offset AGN per 4 
deg^, and even SDSS, with an effective PSE size of 1.4”, could find 
up to 1 per 14 deg^. 

The bottom panels of Eigure[^show the total predicted num¬ 
ber of observable spatially-offset AGN in a given survey, versus 
survey area. Results for spatial offsets are given for T/S'T-COSMOS 
(and a hypot hetical JWST survey with the same area), a Euclid 
wide survey { [Laureijs et al.||201l|), a 2000 deg^ WEIRST weak- 
lensing survey ( jSpergel et al.||2013j ), CDE-N, C-COSMOS, PSl- 
MDS, PSl-3pi, SDSS DR7, and the LSST Main Survey. The ran¬ 
dom spin models predict that 5 - 9 spatially offset AGN may be 
found in HST-COSMOS, while none are likely to be identified with 
Chandra in the CDE-N or C-COSMOS surveys. In contrast, fu¬ 
ture wide-area, high-resolution surveys could detect spatially off¬ 
set AGN in large numbers: WEIRST could detect several thousand 
such objects, and Euclid could detect tens of thousands]^ 

An important prediction of our models is that, because offset 
AGN spend most of their time at large separations, seeing-limited 
surveys should be able to resolve spatially-offset AGN. This sug¬ 
gests that offset AGN may be promising targets for ground-based, 
large-area surveys. For random spins, ~ 15-20 offset AGN may be 
found in PSl-MDS, and hundreds of spatially offset AGN could be 
found in SDSS. The predictions are even more optimistic for LSST 


^ The predictions for Euclid assume a point-source sensitivity of 27.3 AB 
mag in the visible band; the predicted numbers are up to a factor of two 
lower for an assumed sensitivity of 24 AB mag in the J band. 
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Figure 12. The cumulative number of offset AGN observable per deg^ is shown as a function of survey volume (in log z), for flux-limited surveys with 
a given resolution and sensitivity. Top, middle, and bottom panels show results for the random-dry, fgas-b, and cold spin models, respectively. The 
left panels show the number of observable spatially offset AGN. As indicated on the plot legends, the thick solid blue curve corresponds to the sensitivity 
of HST-COSMOS in the F814W band, assuming 0.1” resolution, and the dashed cyan curve corresponds to the target JWST sensitivity at 2 p,m, assuming 
0.07” resolution. The green dotted line corresponds to the target sensitivity of the Euclid wide survey in the visible band, assuming 0.2” resolution. The black 
dash-dotted line corresponds to Chandra observations at the sensitivity of the C-COSMOS survey at 2-10 keV, for 0.5” resolution, and the thin solid orange 
line denotes the PSl-MDS result for the r-band, assuming a seeing-limited resolution of 1”. (The dotted orange line, visible only in the random spin model 
as distinct from the MDS, corresponds to the PSl-3pi survey). Finally, the dot-dot-dashed red line shows results for the SDSS r-band, for a PSF size of 1.4”. 
In each case, the minimum resolvable spatial offset Ai^proj is assumed to be twice the resolution limit. The right panels show predictions for observable 
velocity offset AGN. Only results for HST (solid line) and SDSS (dot-dot-dashed lines) are shown. Gray lines correspond to a minimum resolvable LOS offset 
of Aulos > 600 km s“^. For SDSS, we also show results for Aulos > 1000 km s“^ (magenta line). The random models predict up to a few spatially- 
offset AGN per deg^ detectable with space-based instruments, and up to one velocity-offset AGN per 10 deg^ with Aulos > 1000 km s“^ observable in 
the SDSS. Even in the hybrid models, where spins are nearly aligned in gas-rich mergers, up to a few spatially-offset AGN per 10 deg^ could be found. 
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Figure 13. The cumulative number of observable spatially offset AGN at 2 ; < 3 is shown for various surveys. The top panels show the random spin models 
(random-high and random-dry, from left to right; a small offset along the x-axis is added for clarity). The middle panels show the hybrid spin models 
(fgas-a and fgas-b) in the same manner. The bottom panels show the aligned spin models (hot, cold, and 5 deg). Left column: For each spin model 
and observational survey, the cumulative number of observable recoiling AGN with resolvable (projected) spatial offsets per square degree, out to 2 ; = 3, 
is shown versus the limiting flux in the relevant band. For the optical and NIR bands, the quantity plotted is uFi, for the effective wavelength of the band, 
while the X-ray fluxes correspond to -F 2 -iokeV- Points correspond to various observational surveys, as indicated in the plot legend and as follows. Blue 
triangles: JWST at 2/im, 0.07” resolution. Cyan circles: HST-COSMOS in the F814W band, 0.1” resolution. Green inverted triangles: WFIRST high-latitude 
survey (J band, 0.2” resolution). Green thin diamonds: Euclid wide survey (visible band, 0.2” resolution). Black X’s: C-COSMOS and CDF-N (2-10 keV, 
0.5” resolution). Orange points: a seeing-limited resolution of 1” is assumed for LSST (Main Survey r-band; orange stars) and Pan-STARRS (r-band, for 
PSl-MDS and PSl-3pi; orange squares). Red diamonds: SDSS DRV, assuming a PSF size of 1.4” for the r-band. In each case, the minimum resolvable spatial 
offset Ai^proj is assumed to be twice the resolution limit. Error bars are derived from time-weighted Poisson statistics for the cumulative number counts. 
Right column: the total number of spatially-offset AGN that may be detectable in a given survey is shown, versus survey area. In the same manner as above, 
the points correspond to HST-COSMOS and a hypothetical JWST survey of the same area (1.7 deg^), a 15,000 deg^ Euclid wide survey, a 2000 deg^ WFIRST 
high-latitude survey, the 448 arcmin^ CDF-N and 1.7 deg^ Chandra-COSMOS surveys, the 80 deg^ PSl-MDS, SDSS DRV, the 18,000 deg^ LSST Main 
Survey, and the 30,000 deg^ PSl-3pi survey. The random spin models predict that > 10^ offset AGN could be found with Pan-STARRS, LSST, and WFIRST, 
and > 10^ with Euclid. The hybrid models predict ~ 10^ - 10^ offset AGN to be detectable in these surveys. In the most extreme aligned model (5deg), 
however, recoiling AGN may never be detected. 
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and PSl-3pi, which could detect up to several thousand spatially- 
offset AGN. 

Allowing for some amount of pre-merger BH spin alignment 
reduces the number of observable offset AGN, but does not elim¬ 
inate the possibility of detecting recoils (Figure middle and 
right panels). For the hybrid spin models, less than one offset AGN 
would be expected in HST-COSMOS. However, 2-3 may be de¬ 
tectable in PS 1-MDS, and large populations of spatially offset AGN 
are still detectable in SDSS (up to ~ 200 in DR7), LSST (up to 
~ 700 in the Main Survey), and WFIRST (> 300). The hybrid 
spin models also predict that the PSl-3pi and Euclid wide surveys 
could each detect in excess of 10^ spatially offset AGN. 

The predictions for the three aligned spin models (right pan¬ 
els of Figures differ by several orders of magnitude. For high- 
resolution observations {HST, JWST, Euclid, and WEIRST), the 
hot spin model produces comparable numbers of spatially-offset 
AGN as the hybrid models. A few tens of offset AGN could be 
found with LSST and PSl-3pi, and > 10^ could be detected with 
Euclid. In the cold spin model, a handful might be detectable with 
LSST and PSl-3pi, but wide-area surveys with space-based instru¬ 
ments (i.e., WEIRST and Euclid) would be required to find a pop¬ 
ulation of offset AGN. In the 5 deg model, the maximum displace¬ 
ments achieved barely exceed the angular resolution of SDSS, and 
offset AGN are unlikely to ever be detected - even Euclid is pre¬ 
dicted to find at most a few objects over the entire sky. 

3.5.3 Spectroscopic surveys 

Because SDSS also has spectroscopic coverage, we can estimate 
the number of detectable velocity-offset AGN as well. The ran¬ 
dom spin models produce ~ 0.1 per deg^ velocity offset AGN 
with Aulos > 1000 km s“^ (Eigurep^ left panel), implying that 
nearly 10^ such objects should exist in the DR7 footprint. 

This result is of particular interest, because to date, the SDSS 
spectroscopic quasar database is the only data set on which large- 
scale systematic searches for offset AGN have been carried out. 
[Bonning et al.| ( [2007] ) searched for quasar spectra at 0.1 < z < 
0.81 in which the broad and Mg II lines have symmetric profiles 
and consistent velocity shifts relative to the narrow emission lines. 
They found a null result for Av > 800 km s Using different 
selection criteria, [Tsalmantza et al.| ( [20TT| ) searched for shifted BLs 
in quasar spectra at 0.1 < z < 1.5 and identified 32 candidates of 
interest. [Eracleous et al.| \2012\ searched for BL offsets of Av > 
1000 km s“^ in H/3 for quasars at z < 0.7 and found 88 such 
velocity-offset quasars. However, these BL offsets could also arise 
from high-velocity outflows, double-peaked emitters, or binary BH 
motion, and none have been confirmed to date as recoiling BHs. 

Eor our random spin models, only ~ 10 - 30% of all velocity- 
offset AGN occur at z < 0.7, such that up to several hundred offset 
AGN are predicted in the relevant redshift range. This is still sub¬ 
stantially more than are present in the DR7 quasar catalog, but we 
stress that these observational results cannot be compared at face 
value to our models. There are several important distinctions be¬ 
tween the offset AGN in our sample and the selection criteria for 
SDSS quasars ( [Schneider et al.|20T0] |, which will cause some offset 
AGN to be excluded from the latter. The initial selection of SDSS 
quasar candidates is based on their colors, which means that low- 
luminosity AGN whose fiux is dominated by the host galaxy are 
excluded. 

The SDSS quasars are also required to have a measurable BL 
with EWHM > 1000 km s“^, such that Type II (narrow-line) AGN 
are excluded as well. However, a BL is also required for detecting 


velocity-offset AGN, meaning that some fraction of velocity-offset 
AGN should be undetectable in any case. 

We have also neglected the effects of obscuration in our mod¬ 
els. Obscuring dust in the nuclear region is a very common feature 
of AGN, required to explain the optical, mid-IR, and X-ray spectral 
features of many objects, and may often have a toroidal geometry 
(e.g., |Urry & Padovani| 19951 jUeda et al.|260^[Stern et al.|2005l l. 

While recent work has focused primarily on the role of dusty tori 
in obscuring low-velocity recoil events ( [Raffai et al.[2016| l, we find 
that dust obscuration is unlikely to be a major factor in the ob¬ 
servability of spatially-offset AGN, at least for kick distributions 
that assume some BH spin misalignment. This is because little of 
the obscuring material beyond the BL region should remain bound 
to the BH, and Eigure illustrates that offset AGN are most of¬ 
ten found at separations larger than the scale of the torus (see also 
[Komossa & Merritt|2008[ l. In contrast, velocity offsets will occur 
essentially at the moment of recoil, and thus may be obscured from 
the observer for the early phase of their lifetime. If obscuring ma¬ 
terial extends beyond the dusty AGN torus to galactic scales, ow¬ 
ing perhaps to a coincident nuclear starburst, some recoiling AGN 
could be obscured for a significant portion of their offset lifetimes. 
Recall, however, that a majority of velocity-offset AGN should also 
have spatial offsets detectable with HST, and up to 20% may have 
spatial offsets detectable with SDSS. Thus, obscuration may have 
at most a moderate effect on the observability of velocity-offset 
AGN. 

The SDSS DR7 quasars are additionally required to be fainter 
than z 15 mag (to avoid saturation of the spectrograph) and more 
luminous than Mi = —22. These magnitude limits, at least, have 
only a modest effect on the offset AGN population at z < 0.7, 
if random spins are assumed. If we approximate these cuts for 
our sample (using the same optical bolometric corrections from 
[Hopkins et al.[ ( [2007[ ) that are used throughout the text), the pre¬ 
dicted number of velocity-offset AGN (Av > 1000) is reduced by 
about 30% for the random-dry model, and by only 5% in the 
random-high model. 

An obvious alternative explanation for the observed rarity of 
large BL offsets in SDSS is that pre-merger spin alignment plays 
at least some role in reducing the number of superkicks with Uk> 
1000 km s“^. The hybrid spin models predict that ~ 70 (f gas-b 
model) to 220 (fgas-a) velocity-offset AGN with Av > 1000 
km s“^ may exist in the SDSS footprint, most of which are found 
at z < 0.7. With the estimated magnitude limits described above, 
these numbers are reduced to ~ 30 and 200, respectively. Thus, 
even if spins are preferentially aligned only in gas-rich mergers, 
our results imply that a population of velocity-offset AGN may be 
missed by the SDSS DR7 selection criteria. 

In the hot aligned-spin model, up to 10 velocity-offset AGN 
(with Aulos > 1000 km s“^) could be found with SDSS, with 
only ~ 2 at z < 0.7. None are predicted to be observable in the 
cold or 5deg models. Thus, follow-up studies of candidate re¬ 
coils in SDSS could place indirect constraints on pre-merger BH 
spins. 

3.6 Host galaxy properties 

To characterize the host galaxies of offset AGN, we first examine 
their stellar mass distribution and its evolution with redshift (Eigure 
(M) . The stellar mass distribution for all BH merger hosts peaks at 
~ 10^° Mq at the highest redshifts, and evolves via hierarchical 
growth toward larger masses at lower z. The distribution is limited 
at both the low- and high-mass end by the resolution of the simu- 
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Figure 14. Top row: the distributions of host stellar masses are shown 
for all merging BHs in Illustris (thick black line), in three redshift bins. 
The host stellar mass is defined as the combined mass of the progenitor 
galaxies in the simulation. The thin solid blue line denotes the subset of 
galaxy merger remnants that host an observable spatially offset AGN, for 
the random-dry spin model and assuming HST-COSMOS sensitivity and 
resolution. The thin red dashed line similarly denotes the hosts of observ¬ 
able velocity-offset AGN (Aulos > 600 km s“^). In both cases, only 
offset AGN with a minimum offset lifetime > 10^ yr are shown. Middle 
row: same as top row, but for the fgas-b spin model. Bottom row: same, 
as top row, but assuming SDSS sensitivity and resolution. The host mass 
distribution is broad but peaks at log (M*/Mq) = 10 - 11 for the random 
and aligned spin models. In the hybrid spin models, higher-mass hosts are 
preferred despite their higher escape speeds, owing to their low gas content, 
and observable recoils are very rare at 2 : > 1. 


lation; specifically, the minimum BH mass and stellar/halo particle 
cuts imposed on each progenitor galaxy exclude many dwarf-dwarf 
mergers and satellite-halo mergers. 

For the random spin models, and assuming HST-COSMOS 
sensitivity and resolution, the offset AGN hosts have a broad mass 
distribution that peaks at log (M^/M©) ~ 10 - 11 (Figuretop 
panel). At z < 1, offset AGN favor lower-mass hosts relative 
to the overall distribution. Nominally, low-mass galaxies should 
more easily produce offset AGN owing to their low central escape 
speeds, and indeed, for the random spin models, nearly 30% of BH 
mergers occurring in galaxies with log (M^/M©) < 10 at low red- 
shifts give rise to an observable spatially-offset AGN. 

Ai z > 1, BH merger hosts have lower M* on average, but 
the spatially-offset AGN hosts do not. At most a few percent of 
low-mass galaxy merger remnants (log (M^/M©) < 10) produce 
observable spatial offsets. In part, this is because low-mass galaxies 
generally host smaller BHs with lower AGN luminosites. These 
low-mass bulges are also increasingly gas-rich at high z, making it 
difficult for BHs to escape the nucleus. 
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Figure 15. The distributions of stellar mass (M*) versus the fraction of 
cold, star-forming gas (/gas,sf . relative to M*) are plotted for different red- 
shift bins. The top row shows the distributions for all merging BHs in our Il¬ 
lustris sample (8800 mergers), while the middle and bottom rows show only 
the 492 mergers that produce observable spatially-offset recoiling AGN in 
the random-dry model, assuming HST-COSMOS resolution and sensi¬ 
tivity. In the middle row, the distributions are weighted by number, and 
in the bottom row, they are weighted by time. The overall population of 
merger hosts evolves from low masses and high gas content to higher-mass, 
quenched galaxies. In contrast, the hosts of observable offset AGN favor 
lower gas fractions at all redshifts. 


Figure shows the distribution of M* versus the star form¬ 
ing gas fraction (/gas,sf) for all BH merger hosts and for spatially- 
offset AGN hosts. The merging galaxy population shows a strong 
“quenching” trend evolving from low-mass, gas-rich galaxies at 
high z to higher-mass, gas-poor galaxies at z ~ 0. Comparing 
with the hosts of spatially offset AGN, we see that they follow the 
same general trend, but that galaxies with low- to moderate gas 
fractions are favored; extremely gas-rich systems rarely produce 
offset AGN. This is even clearer in the bottom panel of Figure p~5] 
which shows the same distribution weighted by offset AGN life¬ 
time. The longest-lived offset AGN occur preferentially in galaxies 

with low /gas,sf. 

The M* — /gas,sf distributions for velocity-offset AGN hosts 
are similar to those for spatially offset AGN. Velocity-offset AGN 
occur more often at z > 1 in low-mass, gas-rich galaxies than do 
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spatial offsets, but these generally have short lifetimes, < 1 Myr. 
The suppression of recoils in high- 2 :, gas-rich galaxies is another 
reason that our results depend only weakly on the assumed flux 
limit, because it is primarily at high redshifts where the minimum 
observable flux is more constraining than the minimum /Edd = 
10 “^ 

In the aligned spin models, offset AGN are much rarer over¬ 
all, and almost none occur in the 5deg model. Recoil velocities are 
generally too small to produce observable offsets in either gas-rich 
galaxies or massive, quenched galaxies. As a result, offset AGN in 
the aligned spin models favor similar host stellar masses as the ran¬ 
dom models (log (M^/M©) ~ 10 - 11). No spatially-offset AGN 
are found at 2 : > 1 in the cold spin model, or at 2 : > 2 in the hot 
model. 

When the degree of BH spin alignment depends on the host 
gas fraction, as in the hybrid spin models, recoils are preferentially 
suppressed in gas-rich hosts by spin alignment in addition to the 
aforementioned dynamical effects. As a result, offset AGN are pref¬ 
erentially found in higher-mass hosts, which have lower gas content 
on average. This is consistent with the above observation that mas¬ 
sive BHs dominate the offset AGN population in the hybrid mod¬ 
els. In the f gas-b model (Figuremiddle panels), the extreme 
alignment of the 5 deg spin model for gas-rich hosts eliminates vir¬ 
tually all observable offsets in hosts with log (M^/M©) < 10.5. 
At 2 : > 1, there are no merger remnants with /gas,sf < 0.1, and 
no observable offset AGN. The effect is less severe in the f gas-a 
model, where spins are less highly aligned, but higher-mass, low- 
/gas galaxies are still favored as offset AGN hosts. 

To summarize, we And that observable offset AGN favor host 
masses of log (M*/M©) = 10 - 11 in the random spin models and 
that higher-mass hosts are favored in the hybrid spin models. Offset 
AGN are more likely to be found in hosts with low to moderate cold 
gas fractions, especially in the hybrid spin models. These trends in¬ 
dicate that the host properties of offset AGN could provide indirect 
information about BH spin alignment. However, because the distri¬ 
butions of M* and /gas,sf are broad, a statistical sample of offset 
AGN would be required to distinguish between spin models. Al¬ 
ternatively, flnding even one example of an offset BH in a very 
massive galaxy would strongly argue for spin misalignment in at 
least some cases. 

Distinguishing between BH spin models based on host stel¬ 
lar masses would also likely require space-based observations. The 
bottom panels of Figure show the M* distribution for all BH 
merger hosts and for offset AGN hosts in the random-dry model, 
as in the top panels, but for SDSS resolution and sensitivity. Rel¬ 
ative to the population detectable with HST-COSMOS, there are 
fewer offset AGN overall, and the preference for low-mass hosts 
largely disappears. Because lower-mass BHs (hosted in lower-mass 
galaxies) have shorter recoiling AGN lifetimes, they are more likely 
to exhaust their fuel supply before reaching an offset resolvable 
with ground-based instruments. Also, many low-mass BHs are too 
faint to be detected in SDSS at 2 : > 0.5. 

We again note that the BH merger sample is incomplete in the 
dwarf mass regime (log (M*/M©) < 10) owing to mass resolution 
limits. However, the BH occupation fraction in such galaxies is not 
empirically constrained in any case. At the massive end, the resolu¬ 
tion limits on minor mergers are unimportant, because these events 
do not produce observable recoils. 

It is important to bear in mind that these galaxy merger rem¬ 
nant models are constructed using the progenitor galaxy properties 
in the simulation snapshot prior to each BH merger; as such, they 
do not account for evolution that may occur between the snapshot 


and the time of merger. In particular, the gas fractions shown in 
Figure can be considered “initial” values, in that some of this 
gas will form stars or be expelled during the galaxy merger, prior 
to the BH merger and recoil. 

3.7 Dependence on model parameters 

3. 7 .1 Accretion disk parameters 

Accretion onto recoiling BHs is assumed to occur via a thin, 
radiatively-efflcient “alpha”-disk in the inner regions. At low Ed¬ 
dington ratios, theory predicts that accretion will become radia- 
tively inefficient, and our accretion disk model will become un¬ 
physical. In order to avoid artiflcially long offset AGN lifetimes, 
a minimum Eddington ratio for the recoiling AGN must be im¬ 
posed. This is especially important at low redshift, where AGN ob¬ 
servations are generally not flux-limited. In our flducial models, 
we adopt a minimum f^dd = 10“^ for recoiling AGN, motivated 
by the distribution of Eddington ratios in observed quasars (e.g., 
|Shen & Kelly||2012| |. This is a conservative choice; massive BHs 
in particular may be detectable at lower Eddington ratios. Here we 
explore the consequences of adopting a lower minimum value of 

/Edd = 10“^. 

The statistics of recoil events and properties of their hosts do 
not change; only the observable lifetime of each recoil event de¬ 
pends on the minimum Eddington ratio. Spatially-offset AGN life¬ 
times can be > 100 Myr in the /Edd,min = 10“^ model, while 
the maximum lifetimes for /Edd, min = 10“^ are ~ 40 Myr. Simi¬ 
larly, velocity-offset AGN have a maximum lifetime of ~ 60 Myr, 
versus a maximum of ~ 30 Myr for /Edd,min = 10“^. (Velocity- 
offset AGN lifetimes are more often limited by deceleration in the 
host potential than by the AGN luminosity.) 

In most cases, setting /Edd,min = 10“^ produces a factor 
of a few more spatially offset AGN per square degree than in the 
flducial model. Because the AGN reach larger separations before 
falling below the minimum /Edd, the biggest increase in number 
of observable recoils occurs for seeing-limited observations, where 
spatial resolution is most often a limiting factor. Particularly in 
the aligned spin models, where large spatial offsets are rare, the 
/Edd,min = 10“^ model produces at least 10 times more spatially- 
offset AGN per deg^. Speciflcally, the hot spin model predicts 1-2 
spatially-offset AGN in HST-COSMOS and > 10^ in the LSST 
Main Survey, while even the cold spin model predicts nearly 100 
spatially-offset AGN to be observable with LSST. Velocity-offset 
AGN number counts for /Edd, min = 10“^ are more similar to 
the flducial model, increasing by less than a factor of two. This is 
because at low redshift, the offset AGN lifetimes are still more of¬ 
ten limited by /Edd,min than by the absolute flux limit, even for 
/Edd,min = 10“^. The results are more sensitive to the assumed 
flux limit than in the flducial models, but only moderately so. 

We also consider how our results depend on the choice of 
/riaf = 0.05, the critical Eddington ratio at which accretion tran¬ 
sitions to a radiatively-inefflcient regime. If this regime is ignored 
entirely, such that a constant radiative efficiency of 0.1 is used at 
all times, the maximum spatially-offset AGN lifetime is ~ 90 Myr. 
Similar to the /Edd,min = 10“^ model, these longer lifetimes in¬ 
crease the number of predicted observable recoils by up to a factor 
of a few for most models, and by more than an order of magnitude 
for seeing-limited observations in the aligned spin models. 

Note that assumption of a time-dependent accretion rate is 
especially important in these models, owing to their longer AGN 
lifetimes. Eor both the /Edd,min = 10“^ and radiatively-efflcient 


© 0000 RAS, MNRAS 000, 000-000 




Detection of recoiling BHs 21 


random-high, random-dry 


10 " 


bJO 

q; 


10 '^ 


CO 

V 

10 '^ 


a 


10 

10 " 


-3 


▲ 


t 

\ 




A JWST X Chandra 

• HST ★ LSST 

T WFIRST ■ Pan-STARRS 

♦ Euclid ♦ SDSS 


10 


17 


10 '^® 10 

■^band,lim 

fgas-a, fgas-b 


15 


[erg s ^ cm 


10 ' 

-2l 



hot, cold, 5deg 



^band.lim [erg S ^ CHI 


Figure 16. Same as the left column of Figure but for an a viscosity 
parameter of 0.03. The predicted rates of observable spatially-offset AGN 
are generally a factor of a few higher than in the fiducial model, but in the 
aligned spin models the predictions for seeing-limited observations increase 
by more than an order of magnitude. Similar enhancements relative to the 
fiducial model are seen for the /Edd,min = 10“^ and radiatively-efficient 
models. 


models, assuming a constant accretion rate with AGN lifetime 
Mdisk/M bh reduces the predicted number of spatially-offset AGN 
by a factor of a few to a few tens in most cases. Velocity offset AGN 
lifetimes are also generally shorter by a factor of a few if a constant 
accretion rate is assumed. 

Finally, the a viscosity parameter is also important in deter¬ 
mining the accretion disk lifetime. Our fiducial model assumes 
a = 0.3, motivated by models of self-gravitating accretion disks 
(e.g., |Goodman||2003| l. The nature of the viscous accretion fiow 


is uncertain, however, and a is commonly assumed to be lower. 
Accordingly, we have also tested a model with a = 0.03. Figure 
shows the resulting predictions for the number of observable 
spatially-offset AGN per deg^, in the same manner as Figure 
The lower value of a yields longer offset AGN lifetimes and higher 
average luminosities. The predicted number of observable recoils 
is higher by a factor of a few to > 10 than in the fiducial model, a 
slightly larger enhancement than seen in the /Edd,min = 10“^ or 
radiatively-efficient models. Again, predictions for seeing limited 
observations in the aligned spin models are most sensitive to the 
change in accretion disk parameters. 

Overall, we find that our results depend on the accretion model 
for recoiling AGN, but because the fiducial model is designed to be 
conservative, variations in the accretion model parameters tend to 
increase the observability of recoils. Moreover, the aligned spin 
models are the most sensitive to changes in the accretion rate; 
for the random and hybrid models, the alternate accretion models 
tested here generally increase the predicted number of observable 
offset AGN by factors of a few at most. 


3.7.2 Circumnuclear gas disk 

Our fiducial galaxy models include a compact circumnuclear gas 
disk with a mass equal to the total mass in cold, star-forming gas in 
the progenitor galaxies. This is well-motivated by the correlation 
of galaxy mergers with both nuclear starbursts and AGN fueling, 
which is observed empirically and in simulations. Simulations have 
also shown that the steep central potentials created during gas-rich 
mergers can greatly suppress recoil trajectories ( [Blecha et al.|2011[ 
ISijacki et al.|2011| |. As discussed in Section [T^ the gas component 
included in our galaxy models does indeed reduce offset AGN life¬ 
times in gas-rich galaxies, such that recoils are more likely to be 
detected in merger remnants with relatively low gas fractions. 

However, the simplified nature of these potential models can¬ 
not fully capture the complexity of a galaxy merger remnant. In 
particular, while recoiling BH orbits will generally be centrophilic, 
irregularities in the potential will prevent them from returning ex¬ 
actly to the galaxy center, thus delaying their orbital decay. We have 
also neglected the consumption of cold gas by star formation in the 
interval between the simulation snapshot time (when the progen¬ 
itor galaxy properties are obtained) and the BH merger time, as 
well as during the recoil event. The amount of cold gas present in 
the merger remnant model may therefore be systematically over¬ 
estimated. In order to understand how these assumptions regard¬ 
ing the gas disk model infiuence our results, we consider the case 
in which the gas disk is neglected entirely in the galaxy potential 
model. (The host gas fraction is still used to determine pre-merger 
BH spin alignment in the hybrid spin models.) 

Naturally, galactic escape speeds are lower in the absence of 
a gas disk component. For the random spin models, the escape 
fractions are significant: 20% and 30% of recoil events for the 
random-dry and random-high models, respectively. More¬ 
over, even when scaling the kick speed to the escape speed, return 
times (and maximum displacements) are larger for a given Uk/Uesc 
in the model without gas. Recoils with Uk/Uesc ~ 0.7 can displace 
the BH from the galactic center for up to 1 Gyr, a factor of a few 
longer than in the fiducial model (Figure]^. Recoils with low to 
moderate Uk/Uesc (^ 0.6) are most strongly affected by the shape 
of the inner density profile. When a dense gas component is in¬ 
cluded, some BHs with Uk/Uesc< 0.6 have return times < 10^ yr 
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Figure 17. Same as Figurebut for a model in which no gas disk is in¬ 
cluded in the merger remnant potential. The 5deg spin model is shown. 
Recoil events with Uk/uesc ^ 0.6, which dominate the offset AGN popula¬ 
tion in this model, are most affected by the absence of a gas disk component. 


(Figure [^, but in the model without a gas disk, such events are rare 
even in the aligned spin models. 

As a result, offset AGN lifetimes are longer in the model with 
no gas disk, and a larger fraction of recoil events produce resolv¬ 
able spatial offsets. Figure p^sh ows the same distributions of offset 
AGN properties as in Figure] ll) for the 5deg spin model and HST 
resolution and sensitivity, but for the model with no gas disk. The 
longer return times for moderate-velocity recoils create a popula¬ 
tion of flux-limited offset AGN at Ai^proj ^0-5” that is largely 
absent in the flducial model. These offset AGN would be resolv¬ 
able with HST, JWST, Euclid, and WFIRST; the no-disk model pre¬ 
dicts that up to 14 spatially-offset AGN per deg^ could be found in 
HST-COSMOS, and 7 per deg^ with Euclid and WEIRST. Even the 
cold spin model predicts that 1-2 spatially-offset AGN could be de¬ 
tected in HST-COSMOS. Eor coarser angular resolutions (> 0.5”), 
though, the predicted numbers of spatially-offset AGN are closer 
to the flducial model predictions, because these offsets mostly re¬ 
sult from larger kicks that are less strongly affected by the inner 
potential. 

The offset AGN that are observable in the no-disk model but 
not the flducial model are almost exclusively hosted in gas-rich 
mergers (/gas,sf > 0.3), and most occur at high redshift (z > 1). 
Because gas-rich galaxies have lower masses, on average, low-mass 


hosts are more heavily favored in the no-disk model. Also, in con¬ 
trast to the flducial model, there is no strong preference in the no¬ 
disk model for offset AGN to be hosted in gas-poor galaxies. 


3.7.3 Gas-rich merger definition 

We distinguish between gas-rich and gas-poor galaxy mergers to 
determine which BH binaries are likely to have an ample reservoir 
of cold gas capable of aligning the BH spins. Gas-poor mergers are 
also assumed to require a longer timescale for the BHs to merge, 
though this has little effect on our results relative to the mass-ratio 
dependence of our merger-timescale prescription. We have adopted 
a conservative deflnition of gas-rich mergers as those with a star¬ 
forming gas fraction /gas,sf > 0.1. According to this deflnition, 
more than 85% of mergers at < 1, and nearly all of those at 
z > 1, are gas-rich. In the hybrid spin models, a less stringent 
deflnition of gas-rich mergers would result in fewer aligned spins 
and a greater number of superkicks. 

We test this assertion by considering a model with a critical 
gas fraction of 0.2, for which ~ 75% of galaxies at z < 1 are de- 
flned as gas-rich. Eor the random and aligned spin models, where 
the critical gas fraction enters only into the merger delay timescale, 
the results do not change signiflcantly. The predicted numbers of 
observable recoils differ by a few percent from the flducial model. 
Eor the hybrid spin models, the model with critical /gas,sf = 0.2 
predicts a factor of 2 - 3 more observable recoiling AGN (for spa¬ 
tial and velocity offsets) than the flducial model. The preference 
for high-mass host galaxies also largely disappears, because greater 
numbers of low-mass galaxies are deflned as gas-poor, and thus can 
produce high-velocity kicks. 

In reality, BH spin alignment should depend on the mass (and 
angular momentum) of gas in a circumnuclear disk around the BHs, 
rather than the total gas mass in the galaxy. However, because it is 
unclear what fraction of the cold gas will condense into a coherent 
disk during the merger, we have opted to tie spin alignment to the 
global character of the merger (gas-rich or gas-poor). Eor complete¬ 
ness, we have also tested a model in which the total star-forming 
gas mass is compared to the total BH mass, with gas-rich galax¬ 
ies deflned as those with Mgas,sf/MBH > 10 (such that Mdisk 
> Mbh if ^ 10% of the total cold gas mass is condensed into a 
circumnuclear disk). In this case, > 95% of mergers are gas rich at 
z < 1, indicating that this is an even more conservative deflnition 
of gas-rich mergers than our flducial model. 

Again, for random and aligned spins, there is negligible differ¬ 
ence between this model and the flducial case. Eor hybrid spin mod¬ 
els, despite the small fraction of mergers that do not have aligned 
BH spins, only a factor of ~ 2 fewer offset AGN are predicted with 
this deflnition of gas-rich mergers. We also note that, while offset 
AGN tend to have high BH masses in the hybrid spin models re¬ 
gardless, this effect is even stronger when spin alignment is tied to 
the BH mass. 


3.7.4 BH merger timescales 

As discussed in Section[2?T] a delay is added to the time of each BH 
merger in our model, because BH binary evolution cannot be mod¬ 
eled on sub-grid scales, and because minor-merger timescales in 
particular may be underestimated by the BH prescriptions in Illus- 
tris. If we examine the results without this delay timescale added, 
we And relatively minor differences; the delay timescale lowers the 
predicted numbers of offset AGN by less than a factor of two in 
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most cases. In part, this owes to the small fraction of major merg¬ 
ers for which the BH binary is still unmerged at z = 0 according 
to this prescription (2.5% of mergers with log q > —1.5). Also, 
the unmerged binaries that are removed from our sample are bal¬ 
anced to some extent by higher-redshift binaries that have long de¬ 
lay timescales, such that the corresponding recoil events occur at 
lower redshift, where they are more easily observable. 

As shown in Figure the delay timescale shifts the peak 
of BH merger activity to slightly lower redshift, from z ^ 2 to 
z ~ 1.5. Although this does not produce a large effect on the ob¬ 
servability of recoils, these longer merging timescales may have a 
more significant impact on GW signals from BH mergers. The in¬ 
spiral timescales of BH binaries on subgrid scales, and the resulting 
GW signal, will be explored in detail in Kelley et al. (in prepara¬ 
tion). 


4 DISCUSSION 

4.1 Detection and confirmation of recoiling AGN 

Our models predict the numbers of spatially or kinematically off¬ 
set AGN that are a) resolvable and b) above the sensitivity limit 
for a given observational survey. Once a candidate recoiling AGN 
is identified via a spatial or kinematic offset, however, extensive 
multi-wavelength follow-up observations may be required to rule 
out alternate scenarios. Alternative explanations for existing recoil 
candidates include dual or binary BHs in which only one BH is ac- 
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explanations for the velocity-offset quasar J0927-I-2943 are consis¬ 
tent with recent observations. 


Type Iln SNe may resemble AGN for months or years after 
the explosion; they commonly have both broad and narrow emis¬ 
sion lines and slow decay times owing to interactions of shocked 
ejecta with circumstellar material. However, as their spectral fea¬ 
tures evolve and they eventually fade, such events should be dis¬ 
tinguishable from AGN via continued monitoring and multiwave¬ 
length observations. 

Extreme double-peaked emitters in which only one peak is 
visible may also resemble recoiling AGN, in that they can produce 
systemic velocity shifts in broad emission lines ( |Chen & Halpern| 
|1989[ |Eracleous et al.|[T9^ . Different BLs often have very dif¬ 
ferent velocity shifts in such cases, while to first order, recoiling 
AGN should exhibit a consistent shift in all broad lines (owing to 
the bulk motion of the AGN). Similarly, the orbital motion of sub¬ 
parsec binary BHs in which only one BH is active may give rise 
to velocity-shifted BLs. In this case, the velocity shift will evolve 
on the orbital timescale, while recoiling AGN should exhibit a con¬ 
stant velocity shift. 

One of the more challenging alternative scenarios to rule out 
for candidate recoils is that of a well-separated (~ kpc scale) BH 
pair in which only one BH is actively accreting. This is the case for 
CID-42, which has two distinct nuclei separated by 2.5 kpc, one 
of which is a confirmed AGN with a BL offset of 1300 km s“^ 
( |Civano et al.|2010||2012[|Novak et al.|2015| l. Note that the mag¬ 
nitude of this spatial and velocity offset is quite consistent with the 


predictions of our models, if spins are not highly aligned (Figure 
1^. The offset BL and lack of evidence for an AGN in the second 
nucleus support the recoil scenario, and indirect arguments sug¬ 
gest that the host galaxy properties are more consistent with a post¬ 
merger phase, but it is difficult to rule out the possibility that a qui¬ 
escent BH is harbored in the second nucleus ( |Blecha et al.|2013| ). 

An additional challenge in confirming spatially-offset AGN is 
that recoil events occur, by definition, in recently-merged galaxies. 
The host may therefore be highly disturbed, such that identifying 
the galaxy centroid is nontrivial. In cases where the BH binary in¬ 
spiral time may be long, for example in gas-poor mergers, the host 
may have a more relaxed morphology by the time the BH merger 
and recoil occurs. However, we have demonstrated that if spins are 
not highly aligned, offsets larger than a kpc should be common. 
In such cases, it is not necessary to pinpoint the host centroid to 
great accuracy in order to determine that the AGN is offset. This 
is crucial for any programme to search for offset AGN in ground- 
based, all-sky surveys, where galactic morphological features often 
will not be resolved in great detail. Similarly, because observable 
recoiling AGN also tend to have large velocity offsets (> 500 - 
1000 km s“^), any ambiguity in the rest-frame redshift of the dis¬ 
turbed host (owing to residual motion of the NL gas, for example) 
is unlikely to preclude the detection of such offsets. 

As noted earlier, our models do not account for the possibil¬ 
ity that a recoiling BH may encounter fresh fuel as it settles back 
to the galactic center, which would produce recoiling AGN with 
much smaller spatial offsets. This could be especially relevant for 
recoils in dry merger remnants, where the BH may undergo long- 
lived, small-amplitude oscillations in a core potential ( [Gualandris] 
|& Merritt|2008[[Lena et al.|2014| l. Regardless, our results suggest 
that if BH spins are not highly aligned, a population of spatially- 
offset recoiling AGN may be resolvable with seeing-limited obser¬ 
vations. 

The most unambiguous recoiling AGN candidates will likely 
be those with both spatial and velocity offsets, as in CID-42, or 
with a projected spatial offset large enough that the AGN is well- 
separated from the stellar light of the host galaxy. Our results indi¬ 
cate that a majority of velocity-offset AGN should also have spa¬ 
tial offsets resolvable with HST, and 5 - 20% of spatially-offset 
AGN resolvable with ground-based observations should also have 
LOS velocity offsets > 600 km s“^. Additionally, up to 50% 
of spatially-offset AGN resolvable with the ground-based imaging 
(and up to 25% of those resolvable with HST) should appear at radii 
beyond the host stellar bulge. Thus, the possibility exists of finding 
hundreds of such objects with Pan-STARRS and LSST, if spins are 
not strongly aligned. 

Recoiling AGN with large spatial offsets also tend to have low 
Eddington ratios, as their accretion rate declines over time. We as¬ 
sume that accretion transitions to a radiatively-inefficient regime 
below /Edd = 0.05; this regime may often be associated with radio 
emission in AGN. Thus, high radio-to-optical fiux ratios in offset 
AGN, or even small-scale radio jets, could be another good indica¬ 
tion of a past recoil event. 

We note that the hybrid and aligned spin models predict few 
recoils to be detectable at z > 1. For the random spin models, this 
is not the case; velocity-offset AGN in particular may be found in 
equal or greater numbers at z > 1 as at lower redshift (Figure p^. 
This raises the interesting possibility of detecting recoiling AGN 
near the cosmic peak of quasar and galaxy merging activity, and 
such detections would place additional indirect constraints on pre¬ 
merger BH spins. Evidence of past BH mergers at high redshift 
would also place stronger limits on BH merging timescales than 
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could be obtained in the local Universe. However, if such events 
are observable, they would be more challenging to confirm with 
follow-up observations than their low-redshift counterparts. 

Our models provide only analytic approximations to the prop¬ 
erties and structure of the host galaxies of recoiling AGN, which 
should highly dynamic during major mergers. The total stellar mass 
and star-forming gas content are readily obtained from the progeni¬ 
tors, however, and we can identify trends in these properties of off¬ 
set AGN hosts. Specifically, although the host stellar mass distribu¬ 
tion is broad, stellar masses of log (M^/M©) ~ 10 - 11 are favored 
in the random and aligned spin models. In the hybrid spin models 
(for which spin alignment occurs only in gas-rich mergers), off¬ 
set AGN occur preferentially in higher-mass hosts, log (M^/M©) 
>11. This is a testable prediction, but a statistical sample of offset 
AGN would be required to differentiate between these distinct but 
overlapping distributions. 

We also find that merger remnants with low to moderate star¬ 
forming gas fractions are preferred as offset AGN hosts, as recoil 
trajectories tend to be quickly suppressed in gas-rich mergers. This 
is especially true of the hybrid spin models, where almost no ob¬ 
servable recoils occur in gas-rich merger remnants. We suggest that 
a systematic search for offset AGN might focus initially on bulge- 
dominated or disturbed galaxies with low to moderate star forma¬ 
tion rates. 

Again, we caution that our sample is incomplete at the low- 
mass end owing to resolution limits, such that offset AGN in low- 
mass or dwarf hosts are largely excluded from our analysis. This is 
an important caveat; although the merger rate and BH occupation 
fraction of dwarf galaxies are poorly constrained, the low escape 
speeds of low-mass galaxies should result in resolvable offsets in 
many if not most recoil events, even in galaxies that are gas-rich. 
The recoil candidate SDSS 1133 is one possible example of this 
( |Koss et al.|2014| ). Low-mass merger remnants may therefore host 
a population of spatially-offset AGN in addition to those predicted 
by our models. 

As a final caveat, our results assume that recoiling AGN gen¬ 
erally resemble “normal,” stationary AGN, in that they retain an 
accretion disk and broad line region and have an SED that is well- 
modeled with standard AGN templates. This seems a reasonable 
first-order assumption, given that the region producing most of the 
AGN emission will typically have orbital speeds Uorb ^ and 
will feel the recoil event as a negligible perturbation. Nonetheless, 
the possibility that recoiling AGN appear qualitatively different 
from stationary AGN cannot be discounted, and until recoil events 
can be unambiguously confirmed, this remains an uncertainty for 
any search for such objects. 

The detection of even one recoiling AGN with a large velocity 
offset (> 1000 km s“^) or spatial offset (> 10 kpc), or alterna¬ 
tively the discovery of a massive galaxy with no central BH, would 
strongly indicate that binary BH spins are not efficiently aligned 
in all cases. Our results suggest that up to hundreds or thousands 
of spatially offset AGN may be detectable with Pan-STARRS and 
LSST. While these numbers are dependent on model assumptions, 
as discussed in Section [T7] we have attempted to be conservative 
in constructing the fiducial models. If such a population of recoil¬ 
ing AGN is identified, then further distinctions between BH spin 
and recoil velocity distributions may be possible. If no convincing 
recoil candidates are found in large-area surveys, the simplest ex¬ 
planation may be that spin alignment is effective for most binary 
BHs, or that spin magnitudes are low. Constraints on the efficiency 
of spin alignment could also provide indirect information about the 
nature of the accretion fiow, specifically whether it is likely to be 


coherent on timescales required to align the BH spins, or if it is 
more often characterized by chaotic, randomly-oriented streams. 


4.2 Comparison with previous work 


Our predicted rates of observable offset AGN are substantially 
higher than those of VM08, the only previous study of recoiling 
AGN in a cosmological framework. Here, we identify some key 
differences that contribute to this discrepancy. 

VM08 used semi-analytic merger tree models to determine the 
BH merger rate and host properties, similarly integrating recoil tra¬ 
jectories for each event in a galaxy potential model. Our predic¬ 
tions for the number of observable spatially-offset AGN per deg^ in 
the random-high spin model (Figurecan be compared with 
their results for HST, JWST, and Chandra. For our fiducial mod¬ 
els, these numbers are comparable to the predictions of the VM08 
“halo-only” models, and they are a factor of > 10 - 100 above their 
“halo+bulge” predictions. Given that our fiducial model includes 
not only a halo and compact stellar bulge but also a dense gas disk, 
the recoiling BH return times are much shorter than in a halo-only 
model. Thus, our model predictions are in fact significantly higher 
than those of VM08. 

Neglecting the gas component in our models increases the dis¬ 
crepancy with the VM08 results for a halo+bulge potential, up to a 
factor of several hundred for HST and JWST. If we include only a 
DM halo potential, ~ 15 - 25 times more spatially offset AGN are 
predicted than in our fiducial model. 

The recoil kick velocity distribution is sensitive to the BH 
mass ratio distribution, but the kick formula of [Baker et al.H2008] ) 
used in VM08 depends even more sensitively on q than is predicted 
in more recent work (e.g.,|van Meter et al.| 20 T 0 }|Lousto et al.| 2012 j l. 


Specifically, the Baker et al. ( |2008| ) formula scales with 77 ^, where 
T] = g/(l + is the symmetric mass ratio, w hile the more recent 

^ 2008 t 


Baker et al. 


formulae include an 77 term. As a result, the 
formula produces fewer high-velocity recoils, and the steeper de¬ 
pendence on q amplifies discrepancies in the underlying distribu¬ 
tions. Using this formula in our models lowers our predictions by a 
factor of ~ 2 . 

The ejected accretion disk model of [Loebj ( |2007| ) used in 
VM08 has a steep dependence on BH mass (oc Thus, the 

disk must be capped at Mdisk = Mbh in most cases and depends 
sensitively on the underlying BH mass distribution. Our fiducial 
models instead assume the disk model of [Blecha & Loebj ( |2008| ), 
which accounts for disk self-gravity and yields typical disk masses 
of a few percent of the BH mass. Furthermore, the accretion rate in 
our models is assumed to decline over time as the disk diffuses out¬ 
ward. This increases the predicted number of spatially-offset AGN 
by a factor of a few relative to a constant-M model. 

If we attempt to compare as closely as possible with the results 
of VM08, using the same potential models, recoil kick formula, and 
ejected accretion disk model, we still predict spatially-offset AGN 
to be observable at rates that are tens to hundreds of times higher 
than in any of their models. The discrepancy is larger when compar¬ 
ing to the halo+bulge models than the halo-only case, presumably 
because the bulge mass is directly scaled to the BH mass and thus is 
sensitive to differences in the underlying distribution. Owing to the 
significant differences in numerical techniques and prescriptions 
for BH seeding, growth, and feedback between the semi-analytic 
models and cosmological simulations, it is not surprising that the 
resulting BH merger and recoil properties are qualitatively differ¬ 
ent. The offset AGN population has a nontrivial dependence on the 
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details of each method, but nonetheless, a few important differences 
can be pinpointed. 

First, as noted in Sec |3.1| the BH mass ratio distribution in II- 
lustris is qualitatively different from that shown in VM08, the latter 
having a sharp deficit of mergers with g > 0.3. It is suggested that 
this deficit owes to the minimum mass threshold (Mtot > 10^ M©) 
imposed, indicating that the BH mass distribution differs qualita¬ 
tively from that in Illustris (at least for merging BHs), with lower- 
mass BHs favored in the latter. VM08 also show for comparison a q 
distribution derived from Monte Carlo sampling of their z = 0 BH 
mass function, which is in much better qualitative agreement with 
the Illustris result. We can also compare our results with the semi- 
analytic model results of |B arausse] \20 1 2\ ; they find fewer g > 0.3 
BH mergers than in Illustris as well, but the discrepancy is much 
smaller. If we consider only recoil events with q < 0.3 in our mod¬ 
els, the rates of observable recoils are reduced by a factor of ~ 5 - 
6. 

While it is possible that g ~ 1 mergers are over-represented in 
Illustris, owing to the uncertainties in merger timescales discussed 
in Section im the inspiral times are most likely to be underesti¬ 
mated for minor mergers. We have compensated for this by impos¬ 
ing a minimum BH mass of 10® M© and a merger delay timescale 
that scales with q. 

If the merging BH mass distribution is also different in VM08 
than in Illustris - specifically if lower masses are favored in the 
former, this would contribute to the discrepant results as well. The 
mass distribution is not explicitly given in VM08, but the reason¬ 
able agreement of the integrated merger rate for Mbh > 10®“^ 
M© between Illustris and the other semi-analytic models ( |Sesana| 
|et al.|2004l|Barausse|2012| l suggests that the merging BH mass dis¬ 
tribution in Illustris does not differ too greatly from those studies, 
at least at the high-mass end and integrated over cosmic time. 

We again note that any differences in the mass and mass ratio 
distributions between Illustris and the models of VM08 are ampli¬ 
fied by the sensitive dependence of the [Baker et"^ ( |2008| ) kick 
formula and the |Loeb| ( |2007j l disk model to these parameters. How¬ 
ever, these factors alone do not appear to account for all of the dif¬ 
ference between our comparison model and the results of VM08. 

The cumulative BH merger rate for Mbh >10® M© agrees 
well between Illustris and semi-analytic models, and the rate peaks 
at similar redshift ( jSesana et al.|2004] |Barausse|2012| ). However, 
as noted in Section [TT] there are fewer high-redshift mergers and 
more low-redshift mergers in Illustris, owing to the different BH 
seed formation, growth, and merger prescriptions. Differences in 
the redshift distribution of mergers are likely a significant contrib¬ 
utor to the higher rates of observable recoils predicted in this work. 
Empirical constraints on the BH merger rate are insufficient to dis¬ 
tinguish between these models, but the agreement of the Illustris 
galaxy merger rate with observations is encouraging ( [RodrigueT] 
[Gomez et al.|2015| l, as is the agreement with the observed BH mass 
function and quasar luminosity function { [Sijacki et al.[2015[ l. 


5 SUMMARY 

We have developed models for AGN that are spatially or kinemati¬ 
cally offset from their host galaxies owing to GW recoil kicks fol¬ 
lowing a major merger. Using data from the Illustris cosmological 
simulations, we determine the characteristics of recoiling AGN and 
their host galaxies over cosmic time. Because recoil events may be 
preferentially suppressed in gas-rich mergers, the use of hydrody¬ 
namic simulations is critical for determining which recoiling AGN 


may be observable. Each BH merger in our sample is assigned: a) 
a recoil kick velocity based on the BH mass ratio and assumed spin 
distribution, b) an accretion disk that remains bound to the recoil¬ 
ing BH, and c) a merged host galaxy model based on the progenitor 
galaxy properties. The recoil trajectory is integrated in this galaxy 
potential, and the AGN luminosity is calculated at each timestep. 

This study is the first to consider hydrodynamic effects on 
recoiling BHs in a cosmological framework, and the first to ex¬ 
amine the effect of pre-merger BH spin alignment on the observ¬ 
ability of recoils. Mergers between BHs with aligned or nearly- 
aligned spins produce dramatically lower recoil velocities than do 
randomly-oriented spins, and misaligned spins are required to pro¬ 
duce superkicks that can eject BHs from galaxies entirely. BH spin 
alignment is predicted to occur via interaction with a circumbinary 
gas disk, though the degree of alignment is not well constrained 
and likely depends on the properties of the gas disk. We therefore 
consider three categories of BH spin models: those with randomly- 
oriented spins, those in which BHs always undergo some amount 
of alignment, and “hybrid” models in which spin alignment occurs 
only in gas-rich mergers. 

Our main results can be summarized as follows. 

• Recoiling AGN do not spend most of their time at small 
spatial and velocity offsets, unless pre-merger BH spins are highly 
aligned. The longest-lived offset AGN are produced by recoils near 
the host escape speed. This occurs because marginally-bound and 
escaping BHs experience the least deceleration as they leave the 
host nucleus, and because marginally-bound recoiling BHs spend 
most of their time at large apocenters. Recoil events with Uesc 
are intrinsically rare and have short AGN lifetimes because little 
gas remains bound to the recoiling BH. 

• Spatially-offset AGN are most commonly found at projected 
separations > 0.1” from the host galaxy, with a tail extending 
beyond 10” (if BH spins are not strongly aligned). This indicates 
that even seeing-limited observations can resolve spatially-ojfset 
AGN, making them promising targets for current and future 
wide-area surveys such as SDSS, Pan-STARRS, LSST, Euclid, 
and WEIRST. These spatial offsets often correspond to physical 
projected separations > 1 kpc, such that they may be detectable 
even if the host galaxy is disturbed and its centroid ill-defined. 

• If BH spins are randomly oriented, our models predict 
that nearly 10 spatially-offset AGN may be detectable in HST- 
COSMOS. The Pan-STARRS 1 Medium Deep Survey could find 
~ 15 — 20 offset AGN, and hundreds could be found in SDSS 
DR7. The all-sky surveys LSST and PSl-3pi could detect up 
to several thousand spatially-offset AGN, as could a wide-area 
WEIRST survey. An all-sky survey with Euclid could detect up to 
several tens of thousands of such objects. 

• The random-spin models also predict hundreds of velocity- 
offset AGN with Aulos > 1000 km s“^ in the SDSS DR7 
footprint. At face value, this is in tension with the findings of 
Bonning et al.[ ( [2007[ ), [Tsalmantza et al.[ ( [2011[ ), and [Eracleous] 
et al.[ ( [2012[ ) that large BL offsets are rare among SDSS quasars. 
Many of the offset AGN in our sample may be low-luminosity 
or obscured AGN excluded from the SDSS quasar catalog, but 
this may also be a preliminary indication that pre-merger spin 
alignment is effective in at least some mergers, or that BH spins 
are low. 
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• A population of recoiling AGN may be detected in current 
and future surveys even if some pre-merger BH spin alignment 
occurs. If spin alignment is always very efficient (within ~ 5 ° 
of the orbital angular momentum), recoils are unlikely to ever 
be detected. But if only moderate spin alignment occurs, up to a 
few tens of spatially-offset AGN could be detected with LSST or 
Pan-STARRS, and if spins are only aligned in gas-rich mergers, 
hundreds of recoiling AGN may be observable in such surveys. In 
the latter case, over 2000 recoiling AGN could be detected with 
Euclid. 

• For comparable assumptions regarding the recoil velocity 
distribution and host galaxy properties, we find a much higher 
incidence of observable recoiling AGN than a previous study 
(VM08). We argue that this discrepancy owes primarily to sys¬ 
tematic differences in the underlying redshift, mass and mass ratio 
distributions of merging BHs in Illustris versus the semi-analytic 
models, as well as the use of more physically-motivated models for 
recoiling AGN accretion disks and an updated recoil kick formula 
in this work. 

• A majority of recoiling AGN with resolvable LOS velocity 
offsets may have simultaneous spatial offsets resolvable with HST, 
and ~ 4 — 20% of velocity-offset AGN detectable with SDSS 
may also have spatial offsets resolvable with SDSS. This indicates 
that high-resolution imaging could be key for confirming any 
velocity-offset recoil candidates. 

• A smaller but non-negligible fraction of spatially-offset 
AGN (< 10%) should simultaneously have LOS velocity offsets 
Au > 1000 km s“^. We also find that up to 50% of spatially-offset 
AGN resolvable with seeing-limited observations (and up to 25% 
of those resolvable with HST) should be at projected separations 
larger than the stellar bulge radius, such that they could be distin¬ 
guished from inspiraling dual BHs. AGN exhibiting simultaneous 
spatial and velocity offsets, or that are well separated from the 
host galaxy stellar light, or both, would be among the strongest 
candidate recoils. 

• The preferred stellar masses of offset AGN hosts differ 
depending on the degree of pre-merger BH spin alignment. If 
the spins are always drawn from the same distribution (random 
or aligned), the host stellar mass distribution peaks around log 
(M*/M 0 ) = 10 - 11. For the hybrid spin models, where the degree 
of BH spin alignment depends on the host gas content, high-mass 
hosts are preferred despite their larger escape speeds, because they 
have lower gas fractions on average and are more likely to yield 
superkicks from misaligned BH spins. These models also yield 
very few observable recoils at z > 1 , owing to the gas-rich nature 
of high-redshift galaxies. 

• In all cases, offset AGN are most likely to be found in hosts 
with low to moderate gas fractions, owing to the dynamical sup¬ 
pression of recoil trajectories by dense circumnuclear gas in gas- 
rich systems. The longest-lived offset AGN occur preferentially in 
galaxies with low gas fractions. 

Our quantitative predictions for the observability of recoil¬ 
ing AGN have some dependence on model assumptions, including 
the host galaxy potential model and accretion disk parameters. As 
discussed in Section |3.7| variation with these parameters is often 
substantially larger than the statistical errors shown in Figure 


Because the fiducial model is conservative by design, reasonable 
variations in model parameters tend to increase the predicted num¬ 
ber of observable offset AGN. However, the results also depend on 
the underlying distributions of BH mass, accretion rate, and merger 
rate in Illustris. Owing to uncertainties in the sub-grid BH prescrip¬ 
tions, we cannot exclude the possibility of systematic bias in the 
merging BH population, but the success of Illustris in reproducing 
key observables of the BH population ( [Sijacki et al.|2015| t, and the 
reasonable level of convergence of our results (Appendix 0 are 
encouraging. 

These results strongly motivate a systematic search for offset 
AGN in current and future large-area surveys. Discoveries of GW 
recoils would confirm this prediction of strong-field gravity, and 
they would provide robust evidence of recent BH mergers, thereby 
constraining the timescales for BH binary inspiral. We demonstrate 
that observations of recoiling AGN could also place indirect con¬ 
straints on the efficiency of pre-merger spin alignment, which in 
turn could provide new evidence regarding the nature of the BH 
accretion flow. 
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APPENDIX A: RESOLUTION CONVERGENCE 

The results presented in Section are based on the highest- 
resolution Illustris simulation (Illustris-1). Here we discuss the 
level of convergence between the three Illustris simulations 
(Illustris-1, Illustris-2, and Illustris-3), which differ in mass reso¬ 
lution. The BH merger rate for the three simulations is shown in 
Figure |A1| (top panel), for all BHs and divided by total BH mass 
bins. There is a trend toward convergence for high-mass mergers 
(Mtot > 10® M©) and at low redshift {z < 0.3), but there is a clear 
dependence of merger rate on mass resolution, with higher merger 
rates in the higher-resolution simulations. 

As described in Section |2.3| we consider only mergers for 
which each progenitor has Mbh > 10® M© and is hosted in a 
galaxy with at least 80 stellar and 300 DM particles. Thus, it is 
natural to expect a lower merger rate in the lower-resolution sim¬ 
ulations, where low-mass galaxies are excluded from our analy¬ 
sis. The bottom panel of Figure [aT] compares the merger rate for 
the three Illustris simulations with a consistent mass cut applied 
to the BH host galaxies in each case, rather than a minimum par¬ 
ticle number. Here the Illustris-3 curves are the same as in the top 
panel, corresponding to a minimum stellar (halo) mass of 6.4 x 10^ 
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Figure Al. Top panel: the BH merger rate in Illustris is shown in the same 
manner as in Figure]^ but here the resolution convergence between Illustris- 
1, 2, and 3 is shown. Only mergers for which each progenitor galaxy has at 
least 300 DM and 80 stellar particles are included. As before, the solid 
black lines denote the total merger rate for our BH sample with a de¬ 
lay timescale added to each merger, as described in the text. The orange 
dashed and blue dot-dashed lines show the merger rate separated by total 
BH mass: 6 < log Mbh < 8 and log Mbh > 8? respectively. In all cases, 
thicker lines denote higher resolution simulations. Bottom panel: same as 
top panel, but with a constant minimum mass imposed on progenitor galax¬ 
ies, rather than a minimum particle number. A minimum stellar (halo) mass 
of 6.4 X 10® Mq (1.2 X 10^^ Mq) is assumed, corresponding to 300 DM 
(80 stellar) particles in Illustris-3, such that this curve is identical between 
the two panels. The convergence improves substantially in this case. 


Mq (1.2 X 10^^ M©), and the Illustris-1 and 2 curves assume the 
same minimum mass. The convergence improves substantially in 
this case. This demonstrates that our fiducial sample from Illustris- 
1 has a higher BH merger rate primarily because more low-mass 
galaxies can be resolved. 

FigurejA^illustrates the level of resolution convergence of the 
predicted number of observable spatially-offset AGN in our fiducial 
models. Again, a constant minimum progenitor halo mass is used 
for the comparison. In general, the convergence is quite reasonable; 
in most cases the differences between the simulations are within the 
Poisson errors. Convergence is poorest for the aligned spin models 
(hot and cold) in seeing limited surveys, where recoils are un¬ 
likely to be observed in any case. The higher values for Illustris-3 
in these models reflect its higher BH merger rate when a consistent 
minimum mass cut is applied in all simulations (Figure [ aT] bottom 
panel). 
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Figure A2. As in the top panels of Figure[^ each panel shows the cumulative number of observable spatially-offset AGN per deg^, out to 2 ; = 3, for various 
surveys as a function of limiting flux in the corresponding band. Different panels correspond to different spin models, as indicated on the plot labels. The 
5deg spin model is not shown, as the number counts are nearly all below the plot range. The points are in the same style as FigurefT^ except that here results 
are compared between Illustris-1, 2, and 3. Point size denotes the simulation resolution; the largest points correspond to Illustris-1. As in the bottom panel of 
Figure [AT] a constant minimum DM and stellar mass are assumed for the progenitor halos of merging BHs in each simulation, rather than a constant minimum 
particle number. As a result, the predictions here for Illustris-1 are lower than in Figure [T3] where lower-mass halos are included. 
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